Despite the availability of advanced numerical capabilities and computational power, complex multi-physical, multiscale systems continue to challenge modeling efforts due to their elusive behaviors and nontrivial couplings. High-fidelity simulations have paved the way for simulating such systems, but the associated computational costs prohibit such tools from being used in iterative routines for engineering design and optimization. While numerical strategies for accelerating complex at-scale multi-physical simulations exist, they often concede predictive accuracy and lack rigorous justification. Alternatively, methods for developing rigorous multiscale models (e.g., upscaling techniques) offer a priori modeling error guarantees, but are primarily used in academic settings due to their requirements in time and specialized expertise for handling analytically intractable derivations. In our previous work, we developed an automated upscaling engine, Symbolica, that uses symbolic computation to avoid these overheads and rapidly develop multiscale models via rigorous homogenization theory. In this work, we extend Symbolica’s model development capabilities to accommodate multi-domain structures and demonstrate Symbolica’s generality in model development by applying it to a real-world problem; namely, to analyze thermal runaway in Li-ion battery modules. We obtain homogenized models that capture the sharp spatiotemporal gradients associated with thermal runaway to within the guaranteed modeling error. Upon numerically validating the models, we demonstrate how Symbolica can be used to reduce the level of effort in rigorous multiscale modeling and model implementation for real-world applications. Figure 1