Spatial and spatiotemporal discretizations

Discretizing SPDEs

GaussianMarkovRandomFields.discretize Function

discretize(𝒟::MaternSPDE{D}, discretization::FEMDiscretization)::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.

The SPDE dimension D must match the intrinsic (manifold) dimension of the discretization. For flat meshes this is simply the spatial dimension; for a surface mesh embedded in 3D (e.g. a triangulated sphere), use MaternSPDE{2} — the Laplacian then acts as the Laplace–Beltrami operator on the surface.

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{D, S, G, I, Q, GI, H, CH}

Holds the data needed to discretize an (S)PDE using a Finite Element Method.

The validated user-facing constructor lives in the GaussianMarkovRandomFieldsFEM extension and requires Ferrite to be loaded. See its docstring for details.

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

GaussianMarkovRandomFields.intrinsic_dim Function

intrinsic_dim(f::FEMDiscretization)

Return the intrinsic (manifold) dimension of the discretization, i.e. the reference dimension of its elements. For flat meshes this equals ndim; for meshes embedded in a higher-dimensional space it is smaller, e.g. 2 for a surface mesh in 3D such as a triangulated sphere.

This is the dimension d that enters the statistical properties of an SPDE discretized on the mesh (e.g. α = ν + d/2 and the variance normalization of the Matérn SPDE). The implementation lives in the FEM extension.

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 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 Vec and delegates to the vector method.

source

evaluation_matrix(model::MaternModel)

Return the evaluation matrix mapping FEM DOFs to the stored observation points. Only available when the model was constructed from points (i.e., observation_points is not nothing).

source

evaluation_matrix(model::MaternModel, points::AbstractMatrix)

Return the evaluation matrix mapping FEM DOFs to arbitrary spatial points. Useful for prediction at new locations.

Example

model = MaternModel(train_points; smoothness = 1)
A_pred = evaluation_matrix(model, test_points)

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{RefTriangle, 1}() # hide
qr = QuadratureRule{RefTriangle}(2) # hide
disc = FEMDiscretization(grid, ip, qr)
X = [Vec(0.11, 0.22), 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 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{RefTriangle, 1}() # hide
qr = QuadratureRule{RefTriangle}(2) # hide
disc = FEMDiscretization(grid, ip, qr)
X = [Vec(0.11, 0.22), 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 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,
    cell_coords = nothing,
)

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. Either a constant, or a function of the spatial coordinate returning a velocity vector; in the latter case cell_coords must be passed.

• cell_coords=nothing: The coordinates of the current cell (from getcoordinates), required to evaluate a position-dependent velocity at the quadrature points.

source

GaussianMarkovRandomFields.lump_matrix Function

lump_matrix(A::AbstractMatrix, ::Lagrange{S, 1}) where {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 linear state-space model. Concrete subtypes (and their construction logic) live in the GaussianMarkovRandomFieldsFEM extension.

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

State-space model for the implicit Euler discretization of an SDE. The user-facing constructor lives in the GaussianMarkovRandomFieldsFEM extension.

source

GaussianMarkovRandomFields.ImplicitEulerJointSSMMatrices Type

ImplicitEulerJointSSMMatrices

Joint SSM transition matrices for the implicit Euler scheme.

source