Convective and diffusive operators are discretized such that their symmetries are preserved. The resulting discretization inherits all symmetry-related properties of the continuous formulation. It is shown that a symmetrypreserving discretization is unconditionally stable and conservative. A fourth-order, symmetry-preserving discretization method is developed and tested for the numerical simulation of turbulent (flow and) heat transfer in a channel with surface-mounted cubes, where the temperature is treated as a passive scalar. The Reynolds number (based on the channel width and the mean bulk velocity) is Re=13,000. The results of the numerical simulation agree well with available experimental data.