Understanding condensed-phase polariton experiments requires accurately accounting for both realistic cavity geometries and the interplay between polaritons and material dark modes arising from microscopic molecular interactions. The finite-difference time-domain (FDTD) approach numerically propagates classical Maxwell's equations in the time domain, offering a versatile scheme for modeling polaritons in realistic cavities. However, the simple dielectric functions routinely used in FDTD often fail to describe molecular details. Consequently, standard FDTD calculations, to date, cannot accurately describe processes involving the complex coupling between polaritons and dark modes, such as polariton relaxation, transport, and condensation. For more faithful simulations of the energy flow between the polaritons and dark modes, herein, local bath degrees of freedom coupled to the material polarization are explicitly included in FDTD to describe the dark-mode dynamics. This method -- FDTD with auxiliary bath fields (FDTD-Bath) -- is implemented in the open-source MEEP package by adding a Lorentz-Bath material susceptibility, where explicit bath modes are coupled to conventional Lorentz oscillators. With this Lorentz-Bath susceptibility, linear polariton spectra and Rabi-splitting-dependent polariton relaxation rates in planar Fabry--Pérot cavities are reproduced more accurately than those with the conventional Lorentz susceptibility. Supported by a user-friendly Python interface and efficient MPI parallelism, the FDTD-Bath approach implemented in MEEP is ready to model a wide range of polariton phenomena involving realistic cavity geometries.