  Expanding nebulae are produced by mass loss from stars, especially during
late stages of evolution. Multi-dimensional simulation of these nebulae
requires high resolution near the star and permits resolution that decreases
with distance from the star, ideally with adaptive timesteps. We report the
implementation and testing of static mesh-refinement in the
radiation-magnetohydrodynamics code PION, and document its performance for 2D
and 3D calculations. The bow shock produced by a hot, magnetized, slowly
rotating star as it moves through the magnetized ISM is simulated in 3D,
highlighting differences compared with 2D calculations. Latitude-dependent,
time-varying magnetized winds are modelled and compared with simulations of
ring nebulae around blue supergiants from the literature. A 3D simulation of
the expansion of a fast wind from a Wolf-Rayet star into the slow wind from a
previous red supergiant phase of evolution is presented, with results compared
with results in the literature and analytic theory. Finally the wind-wind
collision from a binary star system is modelled with 3D MHD, and the results
compared with previous 2D hydrodynamic calculations. A python library is
provided for reading and plotting simulation snapshots, and the generation of
synthetic infrared emission maps using torus is also demonstrated. It is shown
that state-of-the-art 3D MHD simulations of wind-driven nebulae can be
performed using PION with reasonable computational resources. The source code
and user documentation is made available for the community under a BSD3
licence.
