How to make a Balance law

Defining the set of solved PDEs in ClimateMachine revolve around defining a BalanceLaw. A balance law solves equations of the form:

∂Y
-- = - ∇ • ( F_{first_order}(Y) + F_{second_order}(Y, Σ) ) + S_{non_conservative}(Y, Σ)
∂t

Here, Y, Σ, F_{first_order}, F_{second_order}, and S_{non_conservative} can be thought of column vectors^[1] expressing:

Y the prognostic state variables, or unknowns of the PDEs to be solved
Σ the gradients of functions of the prognostic state variables
F_{first-order} contains all first-order fluxes (e.g. not functions

of gradients of any variables)

F_{second-order} contains all second-order and higher-order fluxes

(e.g. functions of gradients of any variables)

S_{non_conservative} non-conservative sources^[2]

In order to alleviate users from being concerned with the burden of spatial discretization, users must provide their own implementations of the following methods, which are computed locally at each nodal point:

Variable name specification methods

The vars_state method is used to specify the names of the variables for the following state variable types:

Type	Purpose
`Prognostic`	the variables in the prognostic state vector, typically mass, momentum, and various tracers.
`Auxiliary`	any variables required for the balance law that aren't related to derivatives of the state variables (e.g. spatial coordinates or various integrals) or those needed to solve expensive auxiliary equations (e.g., temperature via a non-linear equation solve)
`Gradient`	the gradients of functions of the prognostic state variables. used to represent values before and after differentiation
`GradientFlux`	the gradient fluxes necessary to impose Neumann boundary conditions. typically the product of a diffusivity tensor with a gradient state variable, potentially equivalent to the second-order flux for a prognostic state variable
`UpwardIntegrals`	any one-dimensional vertical integrals from bottom to top of the domain required for the balance law. used to represent both the integrand and the resulting indefinite integral
`DownwardIntegrals`	any one-dimensional vertical integral from top to bottom of the domain required for the balance law. each variable here must also exist in `vars_state` since the reverse integral kernels use subtraction to reverse the integral instead of performing a new integral. use to represent the value before and after reversing direction

Methods to compute gradients and integrals

Method	Purpose
`compute_gradient_argument!`	specify how to compute the arguments to the gradients. can be functions of prognostic state and auxiliary variables.
`compute_gradient_flux!`	specify how to compute gradient fluxes. can be a functions of the gradient state, the prognostic state, and auxiliary variables.
`integral_load_auxiliary_state!`	specify how to compute integrands. can be functions of the prognostic state and auxiliary variables.
`integral_set_auxiliary_state!`	specify which auxiliary variables are used to store the output of the integrals.
`reverse_integral_load_auxiliary_state!`	specify auxiliary variables need their integrals reversed.
`reverse_integral_set_auxiliary_state!`	specify which auxiliary variables are used to store the output of the reversed integrals.
`update_auxiliary_state!`	perform any updates to the auxiliary variables needed at the beginning of each time-step. Can be used to solve non-linear equations, calculate integrals, and apply filters.
`update_auxiliary_state_gradient!`	same as above, but after computing gradients and gradient fluxes in case these variables are needed during the update.

Methods to compute fluxes and sources

Method	Purpose
`flux_first_order!`	specify `F_{first_order}` for each prognostic state variable. can be functions of the prognostic state and auxiliary variables.
`flux_second_order!`	specify `F_{second_order}` for each prognostic state variable. can be functions of the prognostic state, gradient flux state, and auxiliary variables.
`source!`	specify `S_{non_conservative}` for each prognostic state variable. can be functions of the prognostic state, gradient flux state, and auxiliary variables.

Methods to compute numerical fluxes

Method	Purpose
`wavespeed`	specify how to compute the local wavespeed if using the `RusanovNumericalFlux`.
`boundary_state!`	define exterior nodal values of the prognostic state and gradient flux state used to compute the numerical boundary fluxes.

Methods to set initial conditions

Method	Purpose
`init_state_prognostic!`	provide initial values for the prognostic state as a function of time and space.
`init_state_auxiliary!`	provide initial values for the auxiliary variables as a function of the geometry.

General Remarks

While Y can be thought of a column vector (each row of which corresponds to each state variable and its prognostic equation), the second function argument inside these methods behave as dictionaries, for example:

struct MyModel <: BalanceLaw end

function vars_state(m::MyModel, ::Prognostic, FT)
    @vars begin
        ρ::FT
        T::FT
    end
end

function source!(m::MyModel, source::Vars, args...)
    source.ρ = 1 # adds a source of 1 to RHS of ρ equation
    source.T = 1 # adds a source of 1 to RHS of T equation
end

All equations are marched simultaneously in time.

Reference links

as this can leak conservation of the unknowns and lead to numerical instabilities. It is recommended to use either F_{diffusive} or F_{non_diffusive}, as these fluxes are communicated across elements^[3]

where there may be many [Gauss-Lobatto]^[4] points

1 Column Vectors
2Note that using non-conservative sources should be a final resort,
3MPI communication occurs only across elements, not within each element,
4 Gauss-Lobatto