We present the results of a theoretical investigation of the linear thermal expansion coefficients (TECs) of BeF₂, within a direct Grüneisen formalism where symmetry-preserving deformations are employed. The required physical quantities such as the optimized crystal structures, elastic constants, mode Grüneisen parameters, and phonon density of states are calculated from first-principles. BeF₂ shows an extensive polymorphism at low pressures, and the lowest energy phases [α-cristobalite with space group (SG) P4₁2₁2 and its similar phase with SG P4₃2₁2] are considered in addition to the experimentally observed α-quartz phase. For benchmarking purposes, similar calculations are performed for the rutile phase of ZnF₂, where the volumetric TEC (α_v), derived from the calculated linear TECs along the a (α_a) and c (α_c) directions, is in very good agreement with experimental data and previous theoretical results. For the considered phases of BeF₂, we do not find any negative thermal expansion (NTE). However we observe diverse thermal properties for the distinct phases. The linear TECs are very large, especially α_c of the α-cristobalite phase and its similar phase, leading to giant α_v (∼175 × 10⁻⁶ K⁻¹ at 300 K). The giant α_v arises from large Grüneisen parameters of low-frequency phonon modes, and the C₁₃ elastic constant that is negatively signed and large in magnitude for the α-cristobalite phase. The elastic constants, high-frequency dielectric constants, Born effective charge tensors, and thermal properties of the above phases of BeF₂ are reported for the first time and hence serve as predictions.