  Cosmological simulations require the use of a multiple time-stepping scheme.
Without such a scheme, cosmological simulations would be impossible due to
their high level of dynamic range; over eleven orders of magnitude in density.
Such a large dynamic range leads to a range of over four orders of magnitude in
time-step, which presents a significant load-balancing challenge. In this work,
the extreme adaptivity that cosmological simulations present is tackled in
three main ways through the use of the code SWIFT. First, an adaptive mesh is
used to ensure that only the relevant particles are interacted in a given
time-step. Second, task-based parallelism is used to ensure efficient
load-balancing within a single node, using pthreads and SIMD vectorisation.
Finally, a domain decomposition strategy is presented, using the graph domain
decomposition library METIS, that bisects the work that must be performed by
the simulation between nodes using MPI. These three strategies are shown to
give SWIFT near-perfect weak-scaling characteristics, only losing 25%
performance when scaling from 1 to 4096 cores on a representative problem,
whilst being more than 30x faster than the de-facto standard Gadget-2 code.
