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,
in which 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.