Many low temperature disordered systems are expected to exhibit Poisson-Dirichlet (PD) statistics. In this paper, we focus on the case when the underlying disorder is a logarithmically correlated Gaussian process $φ_N$ on the box $[-N,N]^d\subset\mathbb{Z}^d$. Canonical examples include branching random walk, $*$-scale invariant fields, with the central example being the two dimensional Gaussian free field (GFF), a universal scaling limit of a wide range of statistical mechanics models. The corresponding Gibbs measure obtained by exponentiating $β$ (inverse temperature) times $φ_N$ is a discrete version of the Gaussian multiplicative chaos (GMC) famously constructed by Kahane. In the low temperature or supercritical regime, the GMC is expected to exhibit atomic behavior on suitable renormalization, dictated by the extremal statistics of $φ_N$. Moreover, it is predicted, going back to a conjecture made in 2001 by Carpentier and Le Doussal, that the weights of this atomic GMC has a PD distribution. In a series of works, Biskup and Louidor carried out a comprehensive study of the near maxima of the 2D GFF, and established the conjectured PD behavior throughout the super-critical regime. In another direction, Ding, Roy and Zeitouni established universal behavior of the maximum for a general class of log-correlated Gaussian fields. In this paper we continue this program simply under the assumption of log-correlation and nothing further. We prove that the GMC concentrates on an $O(1)$ neighborhood of the local extrema and the PD prediction holds, in any dimension $d$, throughout the supercritical regime, significantly generalizing past results. Unlike for the 2D GFF, in absence of any Markovian structure for general Gaussian fields, we develop and use as our key input a sharp estimate of the size of level sets, which could have other applications.