jacobian_matrices

Generic utilities to simplify calculation of (contributions to) Jacobian matrices

Represents a term in an evolution equation, which can be composed of multiple sub-terms. The whole right hand side of an evolution equation can be represented as an EquationTerm, as the sum of all the RHS terms.

If the EquationTerm is not a sum, product, or constant then it is a leaf node of the EquationTerm tree that gives a contribution to the Jacobian matrix.

Some slightly clusmy handling of array data (only storing a flattened view of the arrays) is done to ensure that EquationTerm has no type parameters so that the sub_terms::Vector{EquationTerm} has a concrete type.

Note that not all fields can be used for every kind, but all are included (rather than, for example, having multiple subtypes of an abstract EquationTerm type) so that we can make a type-stable tree of EquationTerm objects.

EquationTerm(state_variable::Symbol, array::AbstractArray{mk_float};
             derivatives::Vector{Symbol}=Symbol[],
             upwind_speeds::Vector{T} where T <: Union{Nothing,AbstractArray{mk_float}}=Nothing[],
             second_derivatives::Vector{Symbol}=Symbol[],
             integrand_coordinates::Union{Vector{Tc},Nothing} where Tc <: coordinate=nothing,
             integrand_prefactor::EquationTerm=unit_term,
             dims_coords...)

Construct an EquationTerm object, used to fill Jacobian matrices.

state_variable is the name of the variable in the jacobian_info that this term represents.

array gives the values of this term. nothing may be passed for array, but this is only intended for internal use.

derivatives is a list of any derivatives of state_variable included in this term, which are (optionally) upwinded by upwind_speeds.

second_derivatives is a list of any second derivatives of state_variable included in this term.

integrand_coordinates is a list of coordinates corresponding to dimensions that state_variable is integrated over in this term. May be combined with derivatives. The prefactor of state_variable in the integral is integrand_prefactor.

dims_coords... are name=c::coordinate keyword arguments giving the coordinates corresponding to the dimensions of state_variable.

Jacobian matrix and some associated information.

matrix is a (non-sparse) array containing the Jacobian matrix.

state_vector_entries are Symbols giving the names of the variables in the state vector.

state_vector_numbers is a NamedTuple providing a lookup table for the position of each name in state_vector_entries from the name.

state_vector_dims is a Tuple of Tuples (one for each state vector variable) giving the dimensions of that variable.

state_vector_coords are the mini_coordinate objects corresponding to state_vector_dims.

state_vector_sizes are the total number of points in each state vector variable.

state_vector_dim_sizes are the total numbers of points in each dimension in state_vector_dims.

state_vector_sizes are the total lengths of each variable in the state vector.

state_vector_dim_steps gives the step needed in the flattened state vector to move one point in each dimension in state_vector_dims.

state_vector_offsets gives the offset from indices in the corresponding variable to indices in the full state vector. The first entry of state_vector_offsets is 0. state_vector_offsets has a length one greater than state_vector_entries - the final element is the total size of the state vector.

state_vector_local_ranges gives the index range in each dimension of state_vector_dims that is locally owned by this process - used for parallelised loops over rows (or columns) when dealing with a single variable at a time from the state vector.

row_local_ranges gives an index range of locally owned rows of matrix - used for parallelised loops over rows (or columns) of the full matrix when it is not necessary to know which variable a row corresponds to.

coords is a NamedTuple containing mini_coordinate objects for all the dimensions involved in the jacobian_info.

spectral is a NamedTuple containing mini_discretization_info objects for all the dimensions involved in the jacobian_info.

boundary_skip_funcs is a NamedTuple that gives functions (for each variable in state_vector_entries) that indicate when a grid point should be skipped in the Jacobian because it is set by boundary conditions.

synchronize is the function to be called to synchronize all processes in the shared-memory group that is working with this jacobian_info object.

CompoundTerm(compound_term_expanded::EquationTerm,
             array::AbstractArray{mk_float}; dims_coords...)

Create a 'compound' EquationTerm. Useful when the discrete version of the term is not exactly equal to the combination of its parts, for example a derivative of a product $\partial(ab) = a\partial b + b \partial a$ where the numerical derivative of $ab$ is not exactly equal to the expanded form. The expanded form, passed in as compound_term_expanded, is needed to calculate the functional derivatives (terminology?) that give the Jacobian. The numerical values of this term are passed as array, with the coordinates corresponding to the dimensions of array being the keyword arguments dims_coords....

ConstantTerm(array::AbstractArray; dims_coords...)

Create a constant EquationTerm, which does not depend on any of the state variables. The values of this term are passed as array, with the coordinates corresponding to the dimensions of array being the keyword arguments dims_coords....

NullTerm()

Create a null EquationTerm, to use for example when a term is switched off by some option. Null terms are dropped when summed with non-null terms, or make the whole result null when multiplied by other terms.

add_term_to_Jacobian!(jacobian::jacobian_info, rows_variable::Symbol,
                      prefactor::mk_float, terms::EquationTerm,
                      boundary_speed=nothing)

Add the contribution of terms to jacobian. The terms should be all or part of the evolution equation for rows_variable.

prefactor multiplies all the terms. Usually it will be the timestep.

boundary_speed is the speed needed by jacobian.boundary_skip_funcs to determine whether a grid point is set by the boundary conditions. This will usually be the speed in the z-direction.

add_term_to_Jacobian_row!(jacobian::jacobian_info,
                          jacobian_row::AbstractVector{mk_float},
                          rows_variable::Symbol, prefactor::mk_float,
                          terms::EquationTerm, indices::Tuple,
                          boundary_speed)

Traverse the tree of EquationTerms, accumulating the prefactors until reaching a leaf term (kind = ETsimple) that represents a state vector variable, which makes a contribution to the Jacobian.

create_jacobian_info(coords::NamedTuple, spectral::NamedTuple; comm=comm_block[],
                     synchronize::Union{Function,Nothing}=_block_synchronize,
                     boundary_skip_funcs=nothing, kwargs...)

Create a jacobian_info struct.

kwargs describes the state vector. The keys are the variable names, and the arguments are 2 element Tuples whose first element is the 'region type' for parallel loops over that variable (or nothing if comm = nothing for serial operation) and whose second element Vector or Tuple of Symbols giving the dimensions of the variable.

coords is a NamedTuple giving all the needed coordinates, spectral is a NamedTuple with all the needed 'spectral' objects.

comm is the communicator to use to allocate shared memory arrays. synchronize is the function used to synchronize between shared-memory operations.

boundary_skip_funcs is a NamedTuple whose keys are the variable names and whose values are functions to use to skip points that would be set by boundary conditions (or nothing if no function is needed for the variable).

jacobian_initialize_bc_diagonal!(jacobian::jacobian_info)

Initialize jacobian.matrix to zero, but with ones on the diagonal for rows that correspond to boundary condition points (i.e. those skipped by jacobian.boundary_skip_funcs).

jacobian_initialize_identity!(jacobian::jacobian_info)

Initialize jacobian.matrix with the identity.

jacobian_initialize_identity!(jacobian::jacobian_info)

Initialize jacobian.matrix to zero.