We consider the mean-field analog of the p-star model for homogeneous random networks, and we compare its behavior with that of the p-star model and its classical mean-field approximation in the thermodynamic regime. We show that the partition function of the mean-field model satisfies a sequence of partial differential equations known as the heat hierarchy, and the models connectance is obtained as a solution of a hierarchy of nonlinear viscous PDEs. In the thermodynamic limit, the leading-order solution develops singularities in the space of parameters that evolve as classical shocks regularized by a viscous term. Shocks are associated with phase transitions and stable states are automatically selected consistently with the Maxwell construction. The case p=3 is studied in detail. Monte Carlo simulations show an excellent agreement between the p-star model and its mean-field analog at the macroscopic level, although significant discrepancies arise when local features are compared.