Abstract We study the evolution of a classical harmonic chain with nearest-neighbor interactions starting from domain wall initial conditions. The initial state is taken to be either a product of two Gibbs ensembles (GEs) with unequal temperatures on the two halves of the chain or a product of two generalized Gibbs ensembles (GGEs) with different parameters in the two halves. For this system, we construct the Wigner function and demonstrate that its evolution defines the generalized hydrodynamics (GHD) describing the evolution of the conserved quantities. We solve the GHD for both finite and infinite chains and compute the evolution of conserved densities and currents. For a finite chain with fixed boundaries, we show that these quantities relax as ∼ 1 / t to their respective steady-state values given by the final expected GE or GGE state, depending on the initial conditions. Exact expressions for the Lagrange multipliers of the final expected GGE state are obtained in terms of the steady-state densities. In the case of an infinite chain, we find that the conserved densities and currents at any finite time exhibit ballistic scaling while, at infinite time, any finite segment of the system can be described by a current-carrying non-equilibrium steady state (NESS). We compute the scaling functions analytically and show that the relaxation to the NESS occurs as ∼ 1 / t for the densities and as ∼ 1 / t 2 for the currents. We compare the analytical results from hydrodynamics with those from exact microscopic numerics and find excellent agreement.