A square copper plate is modelled by the explicit two-dimensional diffusion scheme…

PHYSICS · Multi-part · Fermi University Physics · Computational Physics · Scientific programming for physicists · Scientific programming for physicists

Scientific Computing for Physicists, Ch Arrays in memory: layout, cost and the data structures of physical problems

Problem

A square copper plate is modelled by the explicit two-dimensional diffusion scheme met in the first course, uijn+1=uijn+α(ui+1,jn+ui−1,jn+ui,j+1n+ui,j−1n−4uijn)u^{n+1}_{ij} = u^n_{ij} + \alpha\left(u^n_{i+1,j}+u^n_{i-1,j}+u^n_{i,j+1}+u^n_{i,j-1}-4u^n_{ij}\right), on a grid of 401×401401\times401 points with spacing Δx=Δy=5×10−4 m\Delta x = \Delta y = 5\times10^{-4}\,\mathrm{m}. The thermal diffusivity of copper is D=1.1×10−4 m2 s−1D = 1.1\times10^{-4}\,\mathrm{m^2\,s^{-1}} and the time step used is Δt=5×10−4 s\Delta t = 5\times10^{-4}\,\mathrm{s}. The field is stored as a NumPy array u of shape (401,401)(401,401) and dtype float64 (8 bytes per element) in default C order, with the last index contiguous. ((a)) Compute the diffusion number α=DΔt/Δx2\alpha = D\Delta t/\Delta x^2 and state, quoting the stability limit for the explicit scheme on a square two-dimensional grid, whether this run is stable. ((b)) State the strides (in bytes) of u, then give the strides of the interior view u[1:-1, 1:-1] used in the update, and say why this slice is a view rather than a copy and why it is not itself C-contiguous as a standalone block. ((c)) The stencil is written inside the time loop as one array expression, u_new[1:-1,1:-1] = u[1:-1,1:-1] + alpha*(u[2:,1:-1] + u[:-2,1:-1] + u[1:-1,2:] + u[1:-1,:-2] - 4*u[1:-1,1:-1]), executed 50005000 times. Rewrite the update so that no new array is allocated inside the loop, and state which cost of running the loop this rewrite leaves unchanged.

Browse the problem library