PDE Operators

Purpose

The PDEDescription consists of PDEOperators characterising some feature of the model (like friction, convection, exterior forces, optimality conditions etc.) and describe the continuous weak form of the PDE. They can be separated roughly into two categories: linear operators and nonlinear operators.

GradientRobustMultiPhysics.PDEOperatorType
mutable struct PDEOperator{T<:Real, APT<:AssemblyPatternType, AT<:AssemblyType} <: AbstractPDEOperator

common structures for all finite element operators that are assembled with GradientRobustMultiPhysics; better look at the AssemblyPatternType and the constructors

source

The following table lists all available operators and available physics-motivated constructors for them. Click on them or scroll down to find out more details.

Main constructorsSpecial constructorsMathematically
LinearForm$(\mathrm{A}(...), \mathrm{FO}(v))$
BilinearForm$(\mathrm{A}(...,\mathrm{FO}_1(u)), \mathrm{FO}_2(v))$ or $(\mathrm{FO}_1(u), \mathrm{A}(\mathrm{FO}_2(v)))$
LaplaceOperator$(\kappa \nabla u, \nabla v)$
ReactionOperator$(\alpha u, v)$
LagrangeMultiplier$(\mathrm{FO}_1(u), v)$ (automatically assembles 2nd transposed block)
ConvectionOperator$(\beta \cdot \nabla u, v)$ (beta is function or registered unknown)
HookStiffnessOperator2D$(\mathbb{C} \epsilon(u), \epsilon(v))$ (also 1D or 3D variants exist)
ConvectionRotationFormOperator$((a \times \nabla) u, v)$ (a is registered unknown, only 2D for now)
NonlinearForm

Legend: $\mathrm{FO}$ are placeholders for Function Operators, and $\mathrm{A}$ stands for a (linear) Action (its role is explained in the main constructors).

Assembly Type

Many PDE operators need a specification that decides on which set of entities of the mesh (e.g. cells, faces, bfaces, edges) a PDEOperator lives and has to be assembled. This can be steered by the AssemblyType of ExtendableGrids. The AssemblyType can be also used as an argument for Finite Element Interpolations. The following AssemblyTypes are available.

AssemblyTypeDescription
AT_NODESinterpolate at vertices of the mesh (only for H1-conforming FEM)
ON_CELLSassemble/interpolate on the cells of the mesh
ON_FACESassemble/interpolate on all faces of the mesh
ON_IFACESassemble/interpolate on the interior faces of the mesh
ON_BFACESassemble/interpolate on the boundary faces of the mesh
ON_EDGES (*)assemble/interpolate on all edges of the mesh (in 3D)
ON_BEDGES (*)assemble/interpolate on the boundary edges of the mesh (in 3D)
Note

(*) = only reasonable in 3D and still experimental, might have some issues

Custom Linear Operators

It is possible to define custom linearforms and bilinearforms by specifying Function Operators an Action that determines how the operators are combined for the dot-product with the test function operator.

GradientRobustMultiPhysics.LinearFormFunction
LinearForm(
    operator::Type{<:??},
    action::AbstractAction;
    name,
    AT,
    regions,
    factor,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}

generates a LinearForm L(v) = (f,operator(v)) from an action

  • operator : operator applied to test function
  • action : action that computes a result to be multiplied with test function operator

Optional arguments:

  • regions : specifies in which regions the operator should assemble, default [0] means all regions
  • name : name for this LinearForm that is used in print messages
  • AT : specifies on which entities of the grid the LinearForm is assembled (default: ON_CELLS)
  • factor : additional factor that is multiplied during assembly
  • store : stores a vector of the discretised LinearForm with the latest assembly result

Details on the action: The action is an Action consisting of a kernel function with interface (result, input, ...) and additional argument information. During assembly input is ignored (only in this constructor for LinearForms). The result computed by the kernel function is multiplied (dot product) with the operator evaluation of the test function.

source
LinearForm(
    operator::Type{<:??},
    f::GradientRobustMultiPhysics.AbstractUserDataType;
    name,
    kwargs...
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}

generates a LinearForm L(v) = (f,operator(v)) from a DataFunction

  • operator : operator applied to test function
  • action : DataFunction, evaluation is multiplied with test function operator

Optional arguments:

  • regions : specifies in which regions the operator should assemble, default [0] means all regions
  • name : name for this LinearForm that is used in print messages
  • AT : specifies on which entities of the grid the LinearForm is assembled (default: ON_CELLS)
  • factor : additional factor that is multiplied during assembly
  • store : stores a vector of the discretised LinearForm with the latest assembly result
source
LinearForm(operator_test::Type{<:??})
LinearForm(
    operator_test::Type{<:??},
    operators_current::Vector{DataType}
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}
LinearForm(
    operator_test::Type{<:??},
    operators_current::Vector{DataType},
    coeff_from::Vector{Int64}
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}
LinearForm(
    operator_test::Type{<:??},
    operators_current::Vector{DataType},
    coeff_from::Vector{Int64},
    action::AbstractAction;
    regions,
    name,
    AT,
    factor,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, LinearForm, ON_CELLS}

Creates a (PDE description level) LinearForm based on:

  • operator_test : operator for the test function (assumes linearity for that part)
  • operators_current : additional operators for other unknowns
  • coefffrom : either PDE unknown ids or block ids for CurrentSolution given to assemblyoperator! that should be used for operators_current
  • action : an Action with kernel of interface (result, input, kwargs) that takes input (= all but last operator evaluations) and computes result to be dot-producted with test function evaluation (if no action is specified, the full input vector is dot-producted with the test function operator evaluation)

Optional arguments:

  • regions: specifies in which regions the operator should assemble, default [0] means all regions
  • name : name for this LinearForm that is used in print messages
  • AT : specifies on which entities of the grid the LinearForm is assembled (default: ON_CELLS)
  • factor : additional factor that is multiplied during assembly
  • store : stores a vector of the LinearForm with the latest assembly result (e.g. when the operators sits in a system block that has to be reassembled in an iterative scheme)

Details on the action: The action is an Action consisting of a kernel function with interface (result, input, ...) and additional argument information. During assembly input will be filled with the operator evaluations of the other unknowns (i.e. operatorscurrent). The result computed by the kernel function is multiplied (dot product) with the operator evaluation of the test function (i.e. operatortest). If no action is given, the assembly tries to multiply the operator evaluations (that would have been given as input) directly.

source
GradientRobustMultiPhysics.BilinearFormFunction
BilinearForm(
    operators_linear::Vector{DataType},
    operators_current::Vector{DataType},
    coeff_from::Vector{Int64},
    action::AbstractAction;
    name,
    AT,
    APT,
    apply_action_to,
    factor,
    regions,
    transposed_assembly,
    also_transposed_block,
    transpose_factor,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

generates a BilinearForm defined by the following arguments:

  • operators_linear : operator for the two linear arguments (usually ansatz and test function)
  • operators_current : additional operators for other unknowns
  • coefffrom : either PDE unknown ids or block ids for CurrentSolution given to assemblyoperator! that should be used for operators_current
  • action : tells how to further combine the operatorscurrent+operatoransatz evaluations (=input of action) to a result that is multiplied with the test function operator (if no action is specified, the full input vector is dot-producted with the test function operator evaluation)

Optional arguments:

  • applyactionto : specifies which of the two linear arguments is part of the action input ([1] = ansatz, [2] = test)
  • regions : specifies in which regions the operator should assemble, default [0] means all regions
  • name : name for this BilinearForm that is used in print messages
  • AT : specifies on which entities of the grid the BilinearForm is assembled (default: ON_CELLS)
  • APT : specifies the subtype of the APT_BilinearForm AssemblyPattern used for assembly (e.g. for lumping (wip))
  • factor : additional factor that is multiplied during assembly
  • transposed_assembly : transposes the resulting assembled matrix (consider true here for non-symmetric operators)
  • alsotransposedblock : when true toggles assembly of transposed system matrix block
  • transpose_factor : factor for transposed block (default = factor)
  • store : stores a matrix of the BilinearForm with the latest assembly result (e.g. when the operators sits in a system block that has to be reassembled in an iterative scheme)

Details on the action: The action is an Action consisting of a kernel function with interface (result, input, ...) and additional argument information. During assembly input will be filled with the operator evaluations of the other unknowns (i.e. operatorcurrent, if specified) and appended to that the operator evaluation of one of the two linear argument (decided by applyaction_to). The result computed by the kernel function is multiplied (dot product) with the operator evaluation of the other linear argument. If no action is given, the assembly tries to multiply the operator evaluations (that would have been given as input) directly.

source
BilinearForm(
    operators_linear::Vector{DataType}
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}
BilinearForm(
    operators_linear::Vector{DataType},
    action::AbstractAction;
    kwargs...
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

same as other constructor but with operators_current =

source

Special Linear Operators

Below you find special constructors of available common linear and bilinear forms.

GradientRobustMultiPhysics.LaplaceOperatorFunction
LaplaceOperator(

) -> GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}
LaplaceOperator(
    κ;
    name,
    AT,
    ∇,
    regions,
    store
) -> Union{Nothing, GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}}

constructor for a bilinearform that describes a(u,v) = κ (∇u,∇v) where kappa is some constant (diffusion) coefficient.

source
GradientRobustMultiPhysics.ReactionOperatorFunction
ReactionOperator(

) -> GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}
ReactionOperator(
    α
) -> Union{Nothing, GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}, GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}}
ReactionOperator(
    α,
    ncomponents;
    name,
    AT,
    id,
    regions,
    store
) -> Union{Nothing, GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}, GradientRobustMultiPhysics.PDEOperator{Float64, SymmetricBilinearForm, ON_CELLS}}

constructor for a bilinearform a(u,v) = (αu,v) or (u,αv) with some coefficient α that can be a number or an AbstractDataFunction.

source
GradientRobustMultiPhysics.LagrangeMultiplierFunction
LagrangeMultiplier(
    operator::Type{<:??};
    name,
    AT,
    action,
    regions,
    store,
    factor
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform that describes a(u,v) = (A(operator(u)), id(v)) and assembles a second transposed block at the block of the transposed PDE coordinates. It is intended to use to render one unknown of the PDE the Lagrange multiplier for another unknown by putting this operator on the coressponding subdiagonal block of the PDE description.

Example: LagrangeMultiplier(Divergence) is used to render the pressure the LagrangeMultiplier for the velocity divergence constraint in the Stokes prototype.

source
GradientRobustMultiPhysics.ConvectionOperatorFunction
ConvectionOperator(
    a_from::Int64,
    a_operator::Type{<:??},
    xdim::Int64,
    ncomponents::Int64;
    name,
    AT,
    a_to,
    factor,
    ansatz_operator,
    test_operator,
    regions,
    newton,
    store,
    transposed_assembly,
    bonus_quadorder
) -> Union{GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}, GradientRobustMultiPhysics.PDEOperator{Float64, NonlinearForm, ON_CELLS}}

constructs a convection term of the form c(a,u,v) = (aoperator(a)*ansatzoperator(u),test_operator(v)) as a BilinearForm (or NonlinearForm, see newton argument)

  • a_from : id of registered unknown to be used in the spot a
  • a_operator : operator applied to a
  • xdim : expected space dimension
  • ncomponents : expected numer of components of a

optional arguments:

  • newton : generates a NonlinearForm instead of a BilinearForm that triggers assembly of Newton terms for c(u,u,v)
  • ato : position of a argument, set ato = 2 to trigger assembly of c(u,a,v)
  • ansatz_operator : operator used in the spot u (default: Gradient)
  • test_operator : operator used in the spot v (default: Identity)
  • factor : additional factor multiplied in assemblxy (default: 1)
  • regions : specifies in which regions the operator should assemble, default [0] means all regions
  • name : name for this operator that is used in print messages
  • AT : specifies on which entities of the grid the operator is assembled (default: ON_CELLS)
  • store : stores a matrix of the operator with the latest assembly result
source
ConvectionOperator(
    β::GradientRobustMultiPhysics.AbstractUserDataType,
    ncomponents::Int64;
    name,
    store,
    AT,
    ansatz_operator,
    test_operator,
    transposed_assembly,
    regions
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform that describes a(u,v) = ((β ⋅ ∇) u,v) with some user-specified DataFunction β. The user also has to specify the number of components (ncomponents) the convection is applied to. The operators for u and v can be changed (if this leads to something reasonable).

source
GradientRobustMultiPhysics.ConvectionRotationFormOperatorFunction
ConvectionRotationFormOperator(
    beta::Int64,
    beta_operator::Type{<:??},
    xdim::Int64,
    ncomponents::Int64;
    name,
    AT,
    factor,
    ansatz_operator,
    test_operator,
    regions
)

constructor for a bilinearform a(u,v) = (beta x curl(u),v) where beta is the id of some unknown vector field of the PDEDescription, u and v are also vector-fields and x is the cross product (so far this is only implemented in 2D)

source
GradientRobustMultiPhysics.HookStiffnessOperator1DFunction
HookStiffnessOperator1D(
    μ;
    name,
    regions,
    ∇,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform a(u,v) = (μ ∇u,∇v) where C is the 1D stiffness tensor for given μ.

source
GradientRobustMultiPhysics.HookStiffnessOperator2DFunction
HookStiffnessOperator2D(
    μ,
    λ;
    name,
    AT,
    regions,
    ϵ,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform a(u,v) = (C ϵ(u), ϵ(v)) where C is the 2D stiffness tensor for isotropic media in Voigt notation, i.e. ℂ ϵ(u) = 2 μ ϵ(u) + λ tr(ϵ(u)) for Lame parameters μ and λ

In Voigt notation ℂ is a 3 x 3 matrix
ℂ = [c11,c12,  0
     c12,c11,  0
       0,  0,c33]

where c33 = μ, c12 = λ and c11 = 2*c33 + c12

Note: ϵ is the symmetric part of the gradient (in Voigt notation)

source
GradientRobustMultiPhysics.HookStiffnessOperator3DFunction
HookStiffnessOperator3D(
    μ,
    λ;
    name,
    AT,
    regions,
    ϵ,
    store
) -> GradientRobustMultiPhysics.PDEOperator{Float64, BilinearForm, ON_CELLS}

constructor for a bilinearform a(u,v) = (C ϵ(u), ϵ(v)) where C is the 3D stiffness tensor for isotropic media in Voigt notation, i.e. ℂ ϵ(u) = 2 μ ϵ(u) + λ tr(ϵ(u)) for Lame parameters μ and λ

In Voigt notation ℂ is a 6 x 6 matrix
ℂ = [c11,c12,c12,  0,  0,  0
     c12,c11,c12,  0,  0,  0
     c12,c12,c11,  0,  0,  0
       0,  0,  0,c44,  0,  0
       0,  0,  0,  0,c44,  0
       0,  0,  0,  0,  0,c44]   

where c44 = μ, c12 = λ and c11 = 2*c44 + c12

Note: ϵ is the symmetric part of the gradient (in Voigt notation)

source

Examples

Below some examples for operators for custom forms are given:

# Example 1 : div-div bilinearform with a factor λ (e.g. for divergence-penalisation)
operator = BilinearForm([Divergence,Divergence]; factor = λ, name = "λ (div(u),div(v))")

# Example 2 : Gradient jump stabilisation with an item-dependent action and a factor s (e.g. for convection stabilisation)
xFaceVolumes::Array{Float64,1} = xgrid[FaceVolumes]
function stabilisation_kernel(result, input, item)
    result .= input 
    result .*= xFaceVolumes[item]^2
end
action = Action(stabilisation_kernel, [2,2]; dependencies = "I", quadorder = 0 )
operator = BilinearForm([Jump(Gradient), Jump(Gradient)], action; AT = ON_IFACES, factor = s, name = "s |F|^2 [∇(u)]⋅[∇(v)]")

Nonlinear Operators

Nonlinear Operators can be used to auotmatically setup the required Newton terms for the fixpoint algorithms. If the user does not define the jacobian for the kernel function, automatic differentation is used to compute them (with optional sparsity detection), see below for details.

GradientRobustMultiPhysics.NonlinearFormFunction
NonlinearForm(
    operator_test::Type{<:??},
    operators_current::Vector{DataType},
    coeff_from::Vector{Int64},
    action_kernel,
    argsizes::Vector{Int64};
    name,
    AT,
    newton,
    sparse_jacobian,
    jacobian,
    dependencies,
    bonus_quadorder,
    store,
    factor,
    regions
) -> GradientRobustMultiPhysics.PDEOperator{Float64, NonlinearForm, ON_CELLS}

generates a NonlinearForm defined by the following arguments:

  • operator_test : operator for the test function
  • operators_current : additional operators for other unknowns
  • coefffrom : either PDE unknown ids or block ids for CurrentSolution given to assemblyoperator! that should be used for operators_current
  • action_kernel : function of interface (result, input, ...) that computes the nonlinear quantity that should be multiplied with the testfunction operator
  • argsizes : dimensions of [result, input] of kernel function

Optional arguments:

  • dependencies : code String for additional dependencies of the kernel/jacobians (substring of "XTI")
  • jacobian : default = "auto" triggers automatic computation of jacobians by ForwardDiff, otherwise user can specify a function of interface (jacobian, input, ...) with matching dimensions and dependencies
  • sparse_jacobian : use sparsity detection and sparse matrixes for local jacobians ?
  • regions : specifies in which regions the operator should assemble, default [0] means all regions
  • name : name for this NonlinearForm that is used in print messages
  • AT : specifies on which entities of the grid the NonlinearForm is assembled (default: ON_CELLS)
  • factor : additional factor that is multiplied during assembly
  • store : stores a matrix of the discretised NonlinearForm with the latest assembly result
  • bonus_quadorder : increases the quadrature order in assembly accordingly (additional to usual quadorder based on used FESpaces)

Some details: Given some operator G(u), the Newton iteration reads DG u_next = DG u - G(u) which is added to the rest of the (linear) operators in the PDEDescription. The local jacobians (= jacobians of the operator kernel) to build DG needed for this are computed by automatic differentation (ForwardDiff). The user can also specify a jacobian kernel function by hand (which may improve assembly times).

For default dependencies both the kernel functions for the operator and its jacobian have to satisfy the interface

function name(result,input,...)

where input is a vector of the operators of the solution and result is either what then is multiplied with operator2 of the testfunction (or the jacobian).

source

Other Operators

There are some more operators that do not fit into the structures above. Also, in the future, the goal is to open up the operator level for exterior code to setup operators that are assembled elsewhere.

GradientRobustMultiPhysics.FVConvectionDiffusionOperatorType
FVConvectionDiffusionOperator(
    beta_from::Int64;
    μ
) -> FVConvectionDiffusionOperator{Float64}

finite-volume convection diffusion operator (for cell-wise P0 rho)

  • div(μ ∇ρ + β ρ)

For μ = 0, the upwind divergence div_upw(β ρ) is generated For μ > 0, TODO

source
GradientRobustMultiPhysics.DiagonalOperatorType
DiagonalOperator() -> DiagonalOperator{Float64}
DiagonalOperator(
    value::Real;
    name,
    onlynz,
    regions
) -> DiagonalOperator

puts value on the diagonal entries of the cell dofs within given regions

if onlyz == true only values that are zero are changed

can only be applied in PDE LHS

source