  We present a nonequilibrium strong-coupling approach to inhomogeneous systems
of ultracold atoms in optical lattices. We demonstrate its application to the
Mott-insulating phase of a two-dimensional Fermi-Hubbard model in the presence
of a trap potential. Since the theory is formulated self-consistently, the
numerical implementation relies on a massively parallel evaluation of the
self-energy and the Green's function at each lattice site, employing thousands
of CPUs. While the computation of the self-energy is straightforward to
parallelize, the evaluation of the Green's function requires the inversion of a
large sparse $10^d\times 10^d$ matrix, with $d > 6$. As a crucial ingredient,
our solution heavily relies on the smallness of the hopping as compared to the
interaction strength and yields a widely scalable realization of a rapidly
converging iterative algorithm which evaluates all elements of the Green's
function. Results are validated by comparing with the homogeneous case via the
local-density approximation. These calculations also show that the
local-density approximation is valid in non-equilibrium setups without mass
transport.
