Based on a temperature-aware, first-principles-informed phase-field framework, which accounts for thermal fluctuations as well as for the temperature dependence of the thermodynamic potentials, we present results of ultra-high-resolution simulations of the ferroelectric domain evolution in polycrystalline lead zirconate titanate (PZT). Leveraging the parallel efficiency of a Fourier spectral homogenization scheme, we model a micron-sized sample with atomic-unit-cell spatial resolution. In addition, we outline a computational technique to automatically identify and track the different types of domain walls from phase-field data, here applied to large polycrystals containing thousands of grains. We exploit this technique to study the role of domain patterning during polarization reversal under applied electric fields. Results indicate that the density of domain walls and the domain width obey the well-known Kittel-Mitsui-Furuichi-Roitburd square-root law. Moreover, analyzing the statistics of the domain pattern formation of 12,000-grain samples reveals correlations with the grain misorientation, which is in agreement with high-energy x-ray diffraction experiments. Furthermore, we study the occurrence of the two predominant types of domain patterns—monodomain and lamellae/twin domain structures—whose emergence is traced back to the grain orientation. We also highlight the impact of the domain patterns in PZT ceramics on the macroscopically observable piezo- and dielectric material properties.