0

I am developing a computational fluid dynamics solver which computes repeated operations on NumPy arrays. Adjacent elements need to be added and subtracted with each other to calculate spatial derivatives on the grid. The easiest way to do this has been to utilize NumPy's np.roll() function i.e.

uw = 0.5*(grid.u + np.roll(grid.u,1,axis=2) * 0.5*(grid.w + np.roll(grid.w,1,axis=0))

This is problematic as NumPy's roll function creates a new array every time it is called. I am fairly new to advanced NumPy usage and don't have experience with strides. I believe that I just need to return a view of an array rather than a copy for the aforementioned code to be run faster.

I already thought about using a SciPy convolution, but that seems like overkill for what I want to do. GitHub Copilot's solution does not help much in terms of efficiency. I took a look at this answer which seems somewhat useful. However, this only uses fancy indexing on a 2D array (and it seems like there is conflicting info on whether fancy indexing creates a copy or view). I want to do this rolling on an arbitrary axis like np.roll() does and with arbitrary shifting.

Michael Lawrence Garcia
  • 9
  • 2
  • consider an array with a row that is `[0,1,2,3]`, `roll(1)` would produce `[3,0,1,2]`. How could that be a `view`? They didn't write `np.roll` to produce a copy just to do things the slow way! – hpaulj Jul 21 '23 at 00:36
  • The source code can be seen via the docs. Much iof t sets up the indexing for 'arbitrary axis'. But it comes down to `result[res_index] = a[arr_index]`, in other words `fancy indexing`. No one slice can do it. The duplicate writes slices to a new `np.zeros` – hpaulj Jul 21 '23 at 03:17

0 Answers0