[pphysch]

The notion that the simple Laplace solver can scale to grids of arbitrary size, without modification, is a bit misleading for practical purposes. Computational performance will tank or zero out if memory hierarchy constraints are not considered. The author does mention high-performance solutions like multigrid. However, even a basic successive overrelaxation algorithm like the one shown can be partitioned and parallelized, and it is a very good programming exercise to implement a partitioning scheme using MPI or even a low-performance messaging library (or even just optimize for cache sizes on a single device, with no network transport).

Like the transition from the elegant, 5-character Laplace equation to the relatively verbose and complex numerical Julia solver, there is an additional and necessary step in making the numerical solution further scalable with present technology. In particular, the notion that the computational boundaries map nicely to the physical boundaries must be thrown out, because now we must respect the layer of "virtual boundaries" between the partitions.