<div dir="ltr">Hello all,<div><br></div><div>TL;DR: Repa appears to be calculating my stencil operation using boxed arithmetic, which is very slow. How do I make it do unboxed arithmetic?</div><div><br></div><div>I am currently working on a project where I need to simulate a 2D environment that exhibits wave behavior.</div><div><br></div><div>I'm accomplishing this by directly discretizing the wave equation and turning it into a cellular automaton.</div><div><br></div><div>The math here is relatively simple. Let's specify the example to the physical system of a rubber sheet, which exhibits wave behavior. To simulate a rubber sheet, you divide it up into grids. Then, you keep track of the velocity and height of each grid.</div><div><br></div><div>At each iteration of the simulation, you update the velocity of each grid based on the force exerted upon it by its neighboring grids. This is simply a linear multiple of the distance between the grid and its neighbors (i.e. how stretched the sheet is). For example, if cell x is at height X, and its neighbors a,b,c,d are at height A,B,C,D, the tension on X is, for some constant L, </div><div><br></div><div>T = L*((A-X)+(B-X)+...) = L*(A+B+C+D-4X).</div><div><br></div><div>You can calculate this efficiently via a matrix convolution, in particular by a convolution with the matrix</div><div><br></div><div>[[ 0  1   0</div><div>   1  -4  1</div><div>   0  1   0 ]]</div><div><br></div><div>You can be more physically accurate by including the diagonals (weighted down a bit). A popular stencil for such purposes is</div><div><br></div><div><div>[[ 1   2     1</div><div>   2   -12  2</div><div>   1   2     1 ]]</div></div><div><br></div><div>Technically the off-diagonal elements should be something like sqrt(2), but this is sufficient for most purposes.</div><div><br></div><div>After calculating the tension matrix, you divide it by the mass of each grid. This is the acceleration on each grid. Multiply it by the length of your time step, and you get the change in velocity of each grid. Add this to the velocity matrix. Then, multiply the velocity matrix by the time step, and that's the change in height. Add this to the height matrix.</div><div><br></div><div>Boom, you've simulated a single time step of a rubber sheet. Keep doing this over and over to get the time-evolution of the rubber sheet.</div><div><br></div><div>My code to do this is as follows:</div><div><br></div><div><a href="https://gist.github.com/wyager/8c468c9d18ad62aff8bc9738aa947ea4">https://gist.github.com/wyager/8c468c9d18ad62aff8bc9738aa947ea4</a><br></div><div><br></div><div>The imports are messy; sorry about that.</div><div><br></div><div>You'll notice I'm using the Dimensional library, which provides statically typed physical dimensions over any numeric value using a newtype wrapper with typelits. Hopefully this isn't slowing things down somehow. </div><div><br></div><div>According to profiling via "ghc -Odph -rtsopts -threaded -fno-liberate-case -funfolding-use-threshold1000 -funfolding-keeness-factor1000  -optlo-O3 -fforce-recomp -prof -fprof-auto -fexternal-interpreter RepaSim.hs" (whew), fully 80% of time and 84% of allocation comes from the stencil operation. </div><div><br></div><div>The simulation starts with a rubber sheet with a smooth bump in the middle, simulates 100ms of time steps, saves a .png file of the sheet, and then repeats n times (n is typed in via standard input).</div><div><br></div><div>How can I make Repa do this stencil unboxed, and are there other things I should be doing to get better performance?</div><div><br></div><div>Thanks,</div><div>Will</div></div>