The density dependence of nuclear symmetry energy is among the most uncertain parts of the Equation of State (EOS) of dense neutron-rich nuclear matter. It is currently poorly known especially at suprasaturation densities partially because of our poor knowledge about isovector nuclear interactions at short distances. Because of its broad impacts on many interesting issues, to pin down the density dependence of nuclear symmetry energy has been a longstanding and shared goal of both astrophysics and nuclear physics. New observational data of neutron stars including their masses, radii, and tidal deformations since GW170817 have helped improve our knowledge about nuclear symmetry energy especially at high densities. Based on various model analyses of these new data by many people in the nuclear astrophysics community, while our brief review might be incomplete and biased unintentionally, we have learned particularly: (1) The slope parameter L of nuclear symmetry energy at saturation density ρ 0 of nuclear matter from 24 new analyses of neutron star observables is about L ≈ 57.7 ± 19 MeV at 68% confidence level consistent with its fiducial value from surveys of over 50 earlier analyses of both terrestrial and astrophysical data within error bars, (2) The curvature K sym of nuclear symmetry energy at ρ 0 from 16 new analyses of neutron star observables is about K sym ≈ −107 ± 88 MeV at 68% confidence level in very good agreement with the systematics of earlier analyses, (3) The magnitude of nuclear symmetry energy at 2ρ 0 , i.e. E sym (2ρ 0 ) ≈ 51 ± 13 MeV at 68% confidence level, has been extracted from 9 new analyses of neutron star observables consistent with results from earlier analyses of heavy-ion reactions and the latest predictions of the state-of-the-art nuclear many-body theories, (4) while the available data from canonical neutron stars do not provide tight constraints on nuclear symmetry energy at densities above about 2ρ 0 , the lower radius boundary R 2.01 = 12.2 km from NICER's very recent observation of PSR J0740+6620 of mass 2.08 ± 0.07 M and radius R = 12.2 − 16.3 km at 68% confidence level sets a tight lower limit for nuclear symmetry energy at densities above 2ρ 0 , (5) Bayesian inferences of nuclear symmetry energy using models encapsulating a first-order hadron-quark phase transition from observables of canonical neutron stars indicate that the phase transition shifts appreciably both the L and K sym to higher values but with larger uncertainties compared to analyses assuming no such phase transition, (6) The high-density behavior of nuclear symmetry energy affects significantly the minimum frequency necessary to rotationally support GW190814's secondary component of mass (2.50-2.67) M as the fastest and most massive pulsar discovered so far. Overall, thanks to the hard work of many people in the astrophysics and nuclear physics community, new data of neutron star observations since the discovery of GW170817 have enriched significantly our knowledge about the symmetry energy of dense neutron-ric...