  In order to understand the role of QCD in the early universe, we perform
hybrid Monte-Carlo simulation of lattice QCD with $N_f=2+1+1$ optimal
domain-wall quarks at the physical point, on the $64^3 \times
(6,8,10,12,16,20,64)$ lattices, each with three lattice spacings. The lattice
spacings and the bare quark masses are determined on the $64^4$ lattices. The
resulting gauge ensembles provide a basis for studying finite temperature QCD
with $N_f=2+1+1 $ domain-wall quarks at the physical point. In this Proceeding,
we present our first result on the topological susceptibility of the QCD
vacuum. The topological charge of each gauge configuration is measured by the
clover charge in the Wilson flow at the same flow time in physical units, and
the topological susceptibility $ \chi_t(a,T) $ is determined for each ensemble
with lattice spacing $a$ and temperature $T$. Using the topological
susceptibility $\chi_t(a,T) $ of 15 gauge ensembles with three lattice spacings
and different temperatures in the range $T \sim 155-516 $~MeV, we extract the
topological susceptibility $\chi_t(T)$ in the continuum limit.
