We present a general analytic approach to spatially resolve the nano-scale lattice distortion field of strained and defected compact crystals with Bragg coherent x-ray diffraction imaging (BCDI). Our approach relies on fitting a differentiable forward model simultaneously to multiple BCDI datasets corresponding to independent Bragg reflections from the same single crystal. It is designed to be faithful to heterogeneities that potentially manifest as phase discontinuities in the coherently diffracted wave, such as lattice dislocations in an imperfect crystal. We retain fidelity to such small features in the reconstruction process through a Fourier transform -based resampling algorithm designed to largely avoid the point spread tendencies of commonly employed interpolation methods. The reconstruction model defined in this manner brings BCDI reconstruction into the scope of explicit optimization driven by automatic differentiation. With results from simulations and experimental diffraction data, we demonstrate significant improvement in the final image quality compared to conventional phase retrieval, enabled by explicitly coupling multiple BCDI datasets into the reconstruction loss function.