  We introduce an approximation technique for nonlinear hyperbolic systems with
sources that is invariant domain preserving. The method is
discretization-independent provided elementary symmetry and skew-symmetry
properties are satisfied by the scheme. The method is formally first-order
accurate in space. A series of higher-order methods is also introduced. When
these methods violate the invariant domain properties, they are corrected by a
limiting technique that we call convex limiting. After limiting, the resulting
methods satisfy all the invariant domain properties that are imposed by the
user (see Theorem~7.24). A key novelty is that the bounds that are enforced on
the solution at each time step are necessarily satisfied by the low-order
approximation.
