Spatial and spatiotemporal discretizations

Discretizing SPDEs

GaussianMarkovRandomFields.discretize Function

discretize(𝒟::MaternSPDE{D}, discretization::FEMDiscretization{D})::AbstractGMRF where {D}

Discretize a Matérn SPDE using a Finite Element Method (FEM) discretization. Computes the stiffness and (lumped) mass matrix, and then forms the precision matrix of the GMRF discretization.

source

discretize(spde::AdvectionDiffusionSPDE, discretization::FEMDiscretization,
ts::AbstractVector{Float64}; colored_noise = false,
streamline_diffusion = false, h = 0.1) where {D}

Discretize an advection-diffusion SPDE using a constant spatial mesh. Streamline diffusion is an optional stabilization scheme for advection-dominated problems, which are known to be unstable. When using streamline diffusion, h may be passed to specify the mesh element size.

source

Spatial discretization: FEM

GaussianMarkovRandomFields.FEMDiscretization Type

FEMDiscretization(
    grid::Ferrite.Grid,
    interpolation::Ferrite.Interpolation,
    quadrature_rule::Ferrite.QuadratureRule,
    fields = ((:u, nothing),),
    boundary_conditions = (),
)

A struct that contains all the information needed to discretize an (S)PDE using the Finite Element Method.

Arguments

• grid::Ferrite.Grid: The grid on which the discretization is defined.

• interpolation::Ferrite.Interpolation: The interpolation scheme, i.e. the type of FEM elements.

• quadrature_rule::Ferrite.QuadratureRule: The quadrature rule.

• fields::Vector{Tuple{Symbol, Union{Nothing, Ferrite.Interpolation}}}: The fields to be discretized. Each tuple contains the field name and the geometric interpolation scheme. If the interpolation scheme is nothing, interpolation is used for geometric interpolation.

• boundary_conditions::Vector{Tuple{Ferrite.BoundaryCondition, Float64}}: The (soft) boundary conditions. Each tuple contains the boundary condition and the noise standard deviation.

source

GaussianMarkovRandomFields.ndim Function

ndim(f::FEMDiscretization)

Return the dimension of space in which the discretization is defined. Typically ndim(f) == 1, 2, or 3.

source

Ferrite.ndofs Method

ndofs(f::FEMDiscretization)

Return the number of degrees of freedom in the discretization.

source

GaussianMarkovRandomFields.evaluation_matrix Function

evaluation_matrix(f::FEMDiscretization, X)

Return the matrix A such that A[i, j] is the value of the j-th basis function at the i-th point in X.

Arguments

• f::FEMDiscretization: The finite element discretization

• X: Evaluation points, either a Vector of Tensors.Vec or an AbstractMatrix where each row is a point

• field: Field name (default: first field in dof_handler)

Matrix format

When X is a matrix, it should be N×D where N is the number of points and D is the spatial dimension.

source

evaluation_matrix(f::FEMDiscretization, X::AbstractMatrix; field = :default)

Convenience method that accepts a matrix where each row is a point. Converts the matrix to a Vector of Tensors.Vec and delegates to the vector method.

source

GaussianMarkovRandomFields.node_selection_matrix Function

node_selection_matrix(f::FEMDiscretization, node_ids)

Return the matrix A such that A[i, j] = 1 if the j-th basis function is associated with the i-th node in node_ids.

source

GaussianMarkovRandomFields.derivative_matrices Function

derivative_matrices(f::FEMDiscretization{D}, X; derivative_idcs = [1])

Return a vector of matrices such that mats[k][i, j] is the derivative of the j-th basis function at X[i], where the partial derivative index is given by derivative_idcs[k].

Examples

We're modelling a 2D function u(x, y) and we want the derivatives with respect to y at two input points.

using Ferrite # hide
grid = generate_grid(Triangle, (20,20)) # hide
ip = Lagrange{2, RefTetrahedron, 1}() # hide
qr = QuadratureRule{2, RefTetrahedron}(2) # hide
disc = FEMDiscretization(grid, ip, qr)
X = [Tensors.Vec(0.11, 0.22), Tensors.Vec(-0.1, 0.4)]

mats = derivative_matrices(disc, X; derivative_idcs=[2])

mats contains a single matrix of size (2, ndofs(disc)) where the i-th row contains the derivative of all basis functions with respect to y at X[i].

source

derivative_matrices(f::FEMDiscretization, X::AbstractMatrix; kwargs...)

Convenience method that accepts a matrix where each row is a point. Converts the matrix to a Vector of Tensors.Vec and delegates to the vector method.

source

GaussianMarkovRandomFields.second_derivative_matrices Function

second_derivative_matrices(f::FEMDiscretization{D}, X; derivative_idcs = [(1,1)])

Return a vector of matrices such that mats[k][i, j] is the second derivative of the j-th basis function at X[i], where the partial derivative index is given by derivative_idcs[k]. Note that the indices refer to the Hessian, i.e. (1, 2) corresponds to ∂²/∂x∂y.

Examples

We're modelling a 2D function u(x, y) and we want to evaluate the Laplacian at two input points.

using Ferrite # hide
grid = generate_grid(Triangle, (20,20)) # hide
ip = Lagrange{2, RefTetrahedron, 1}() # hide
qr = QuadratureRule{2, RefTetrahedron}(2) # hide
disc = FEMDiscretization(grid, ip, qr)
X = [Tensors.Vec(0.11, 0.22), Tensors.Vec(-0.1, 0.4)]

A, B = derivative_matrices(disc, X; derivative_idcs=[(1, 1), (2, 2)])
laplacian = A + B

source

second_derivative_matrices(f::FEMDiscretization, X::AbstractMatrix; kwargs...)

Convenience method that accepts a matrix where each row is a point. Converts the matrix to a Vector of Tensors.Vec and delegates to the vector method.

source

Utilities

GaussianMarkovRandomFields.assemble_mass_matrix Function

assemble_mass_matrix(
    Ce::SparseMatrixCSC,
    cellvalues::CellValues,
    interpolation;
    lumping = true,
)

Assemble the mass matrix Ce for the given cell values.

Arguments

• Ce::SparseMatrixCSC: The mass matrix.

• cellvalues::CellValues: Ferrite cell values.

• interpolation::Interpolation: The interpolation scheme.

• lumping::Bool=true: Whether to lump the mass matrix.

source

GaussianMarkovRandomFields.assemble_diffusion_matrix Function

assemble_diffusion_matrix(
    Ge::SparseMatrixCSC,
    cellvalues::CellValues;
    diffusion_factor = I,
)

Assemble the diffusion matrix Ge for the given cell values.

Arguments

• Ge::SparseMatrixCSC: The diffusion matrix.

• cellvalues::CellValues: Ferrite cell values.

• diffusion_factor=I: The diffusion factor.

source

GaussianMarkovRandomFields.assemble_advection_matrix Function

assemble_advection_matrix(
    Be::SparseMatrixCSC,
    cellvalues::CellValues;
    advection_velocity = 1,
)

Assemble the advection matrix Be for the given cell values.

Arguments

• Be::SparseMatrixCSC: The advection matrix.

• cellvalues::CellValues: Ferrite cell values.

• advection_velocity=1: The advection velocity.

source

GaussianMarkovRandomFields.lump_matrix Function

lump_matrix(A::AbstractMatrix, ::Lagrange{D, S, 1}) where {D, S}

Lump a matrix by summing over the rows.

source

lump_matrix(A::AbstractMatrix, ::Lagrange)

Lump a matrix through HRZ lumping. Fallback for non-linear elements. Row-summing cannot be used for non-linear elements, because it does not ensure positive definiteness.

source

GaussianMarkovRandomFields.assemble_streamline_diffusion_matrix Function

assemble_streamline_diffusion_matrix(
    Ge::SparseMatrixCSC,
    cellvalues::CellValues,
    advection_velocity,
    h,
)

Assemble the streamline diffusion matrix Ge for the given cell values.

Arguments

• Ge::SparseMatrixCSC: The streamline diffusion matrix.

• cellvalues::CellValues: Ferrite cell values.

• advection_velocity: The advection velocity.

• h: The mesh size.

source

GaussianMarkovRandomFields.apply_soft_constraints! Function

apply_soft_constraints!(K, f_rhs, ch, constraint_noise; Q_rhs = nothing, Q_rhs_sqrt = nothing)

Apply soft constraints to the Gaussian relation

$$\mathbf{K}\mathbf{u}\sim \mathcal{N}({\mathbf{f}}_{\text{rhs}},{\mathbf{Q}}_{\text{rhs}}^{-1})$$

Soft means that the constraints are fulfilled up to noise of magnitude specified by constraint_noise.

Modifies K and f_rhs in place. If Q_rhs and Q_rhs_sqrt are provided, they are modified in place as well.

Arguments

• K::SparseMatrixCSC: Stiffness matrix.

• f_rhs::AbstractVector: Right-hand side.

• ch::ConstraintHandler: Constraint handler.

• constraint_noise::Vector{Float64}: Noise for each constraint.

• Q_rhs::Union{Nothing, SparseMatrixCSC}: Covariance matrix for the right-hand side.

• Q_rhs_sqrt::Union{Nothing, SparseMatrixCSC}: Square root of the covariance matrix for the right-hand side.

source

Temporal discretization and state-space models

GaussianMarkovRandomFields.JointSSMMatrices Type

JointSSMMatrices

Abstract type for the matrices defining the transition of a certain linear state-space model of the form

$$G(\Delta t)x_{k+1}\mid xₖ\sim 𝒩(M(\Delta t)xₖ,\Sigma )$$

Fields

• Δt::Real: Time step.

• G::LinearMap: Transition matrix.

• M::LinearMap: Observation matrix.

• Σ⁻¹::LinearMap: Transition precision map.

• Σ⁻¹_sqrt::LinearMap: Square root of the transition precision map.

• constraint_handler: Ferrite constraint handler.

• constraint_noise: Constraint noise.

source

GaussianMarkovRandomFields.joint_ssm Function

joint_ssm(x₀::GMRF, ssm_matrices::Function, ts::AbstractVector)

Form the joint GMRF for the linear state-space model given by

$$G(\Delta tₖ)x_{k+1}\mid xₖ\sim 𝒩(M(\Delta tₖ)xₖ,\Sigma )$$

at time points given by ts (from which the Δtₖ are computed).

source

GaussianMarkovRandomFields.ImplicitEulerSSM Type

ImplicitEulerSSM{X,S,GF,MF,MIF,BF,BIF,TS,C,V}(
    x₀::X,
    G::GF,
    M::MF,
    M⁻¹::MIF,
    β::BF,
    β⁻¹::BIF,
    spatial_noise::S,
    ts::TS,
    constraint_handler::C,
    constraint_noise::V,
)

State-space model for the implicit Euler discretization of a stochastic differential equation.

The state-space model is given by

G(Δt) xₖ₊₁ = M(Δt) xₖ + M(Δt) β(Δt) zₛ

where zₛ is (possibly colored) spatial noise.

source

GaussianMarkovRandomFields.ImplicitEulerJointSSMMatrices Type

ImplicitEulerJointSSMMatrices{T,GM,MM,SM,SQRT,C,V}(
    ssm::ImplicitEulerSSM,
    Δt::Real
)

Construct the joint state-space model matrices for the implicit Euler discretization scheme.

Arguments

• ssm::ImplicitEulerSSM: The implicit Euler state-space model.

• Δt::Real: The time step.

source