Recalling the expectation of an extremely strong primordial magnetic field $H$, we recheck transitions among the phases of chiral symmetry restoration ($\chi SR$), chiral symmetry breaking ($\chi SB$), and pion superfluidity ($\pi SF$)in the QCD epoch of the early Universe. For homogeneous phases in a finite $H$, a sensible scheme is adopted to determine the phase boundaries of $\pi SF$, which is also superconductivity phase itself. In the first part, the QCD phase diagrams are studied in detail within the chiral effective Polyakov-Nambu--Jona-Lasinio model and the transitions involving $\pi SF$ are found to be of first order at relatively small $H$. As expected from the Meissner effect, the regime of $\pi SF$ shrinks with increasing $H$ and completely vanishes beyond a threshold value. In the second part, the bubble dynamics is illuminated for the stronger first-order transition, $\chi SR\rightarrow\pi SF$, in the more convenient Polyakov-quark-meson model. The coupled equations of motion of pion condensate and magnetic field are solved consistently to give the bubble structure. Then, based on bubble collisions, we explore gravitational wave (GW) emission by developing a simple toy model in advance; and the characteristic frequency of the relic GW is estimated to be of the order $0.1-1\,{\rm K}$ or $10^9-10^{10}\,{\rm Hz}$ in our galaxy.