  At the end of inflation, dynamical instability can rapidly deposit the energy
of homogeneous cold inflaton into excitations of other fields. This process,
known as preheating, is rather violent, inhomogeneous and non-linear, and has
to be studied numerically. This paper presents a new code for simulating scalar
field dynamics in expanding universe written for that purpose. Compared to
available alternatives, it significantly improves both the speed and the
accuracy of calculations, and is fully instrumented for 3D visualization. We
reproduce previously published results on preheating in simple chaotic
inflation models, and further investigate non-linear dynamics of the inflaton
decay. Surprisingly, we find that the fields do not want to thermalize quite
the way one would think. Instead of directly reaching equilibrium, the
evolution appears to be stuck in a rather simple but quite inhomogeneous state.
In particular, one-point distribution function of total energy density appears
to be universal among various two-field preheating models, and is exceedingly
well described by a lognormal distribution. It is tempting to attribute this
state to scalar field turbulence.
