We present an analytical and numerical study of the Ising model on a bilayer honeycomb lattice including interlayer frustration and coupling with an external magnetic field. First, we discuss the exact T=0 phase diagram, where we find finite entropy phases for different magnetizations. Then, we study the magnetic properties of the system at finite temperature using complementary analytical techniques (Bethe lattice) and two types of Monte Carlo algorithms (Metropolis and Wang-Landau). We characterize the phase transitions and discuss the phase diagrams. The system presents a rich phenomenology: There are first- and second-order transitions, low-temperature phases with extensive degeneracy, and order-by-disorder state selection.