Hello, The solution that I mentioned yesterday does work in serial. q0 = Function(Vdg).interpolate(Expression("0.0")) q0.dat.data[:] += 0.01*np.random.randn(*q0.dat.shape) However, when I try running it using mpirun on my mac it fails because the declaration of a random variable is not broadcast. : operands could not be broadcast together with shapes (2500,) (2724,) (2500,) ValueError: operands could not be broadcast together with shapes (2500,) (2908,) (2500,) q0.dat.data[:] += 0.01*np.random.randn(*q0.dat.shape) ValueError: operands could not be broadcast together with shapes (2500,) (2924,) (2500,) Does anyone have a suggestion as to how one can do this to maintain the ability to run this in parallel? Thanks again, Francis ------------------ Francis Poulin Associate Professor Department of Applied Mathematics University of Waterloo email: fpoulin@uwaterloo.ca Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637 ________________________________ From: Francis Poulin Sent: Saturday, April 23, 2016 10:39 PM To: firedrake@imperial.ac.uk Subject: RE: [firedrake] how to set up random initial conditions? Hello again, Thanks Colin, that helps a great deal. q0 = Function(Vdg).interpolate(Expression("0.0")) q0.dat.data[:] = 0.01*np.random.randn(*q0.dat.shape) As I mentioned I want to put together an example that shows 2D turbulence, and this will help a great deal. I also plan to put together an example on the barotropic instability of a jet. If any of these are of interest I would be happy to try and put together a demo or something. Francis