Skip to main content

On the edge

One of the more difficult aspects of coding stencils or other neighbourhood-based operations is working out what should happen at the edge of the grid. Parallel array languages such as ZPL have expended significant effort to provide mechanism to make coding edge effects easier.

PM handles edge effects by making use of optional data types - containers that can hold either a value (e.g.: an integer) or null. Optional types have many uses, including the representation of unbounded data structures in a pointerless language. A number of operators are defined for optional types, most usefully |.  x|y returns the value contained in x if that value is not null. Otherwise y is returned. 

Binary communicating operators (x@[neighbourhood] or x@{neighbourhood}) always return optionally-typed values, with off-edge cells set to null. Combined with the | operator, this makes dealing with edge-effects straightforward.

The following code fragment demonstrates how PM handles edge effects in practice. It solves the 2D steady-state heat transfer problem for a rectangular box whose four walls are maintained at temperatures boundary1..boundary4. A simple solution to this problem involves repeatedly replacing the temperature in each grid cell with the average of its four immediate neighbours until no further change occurs (Jacobi Iteration).

 for cell in [1..N,1..M] do  
      temp := 0.0  
      repeat   
           newtemp := (  
               temp@[-1, 0]|boundary1  
              +temp@[ 1, 0]|boundary2  
              +temp@[ 0,-1]|boundary3  
              +temp@[ 0, 1]|boundary4  
           )/4.0  
         diff := abs(newtemp-temp)  
         temp = newtemp  
      until max::diff < tolerance  
 endfor

This algorithm is a popular example used for discussing and benchmarking parallel systems. Example code is available online for ZPL, MPI and Chapel (amongst others.)

Comments

Post a Comment

Popular posts from this blog

‘Tis the season

… to give an end-of-year update of what is happening with the PM language. It is a long time now since I wrote a PM blog. However, behind the scenes there has been a lot of work done on the guts of the compiler with new inference and optimisation passes added and others taken apart, polished, and put back together again. New language features have also been implemented, including methods, improved object lifetime management and, most importantly, a radically simplified approach to coding stencils. The handling of stencils was one aspect of PM-0.4 with which I was least happy. Since these are a key requirement for nearly all numerical modelling, I have long been looking for a way to make their coding as straightforward as possible. However, this has been one of the greatest challenges in terms of implementation, since the amount of code restructuring need to make a stencil operate efficiently (separating out halo computation, overlapping computation and halo exchange, tiling and the in...
Post a Comment
Read more

The PM Type System

In common with other aspects of PM design, the type system is designed to facilitate the use of flexible programming constructs, such as polymorphism and dynamic dispatch, while emphasising the generation of fast, static, code.  In common with many other modern languages, PM also attempts to combine the safety of static typing with the expressivity of dynamically typed languages such as Python . The type of any PM expression may usually be determined by a static examination of the code. Variables take their types from an initialising expression using the ‘:=’ declaration syntax popularised by Go . a:= 100 ! Declare ‘a’ as an integer (int) variable Variables do not change type during the course of their lifetime – polymorphic programming requires the use of special values, described later. Composite values such as structures are generated using specialised expressions rather than through the invocation of type-specific creator functions: b:= struct var_descriptor...
Post a Comment
Read more

Communicating Operators - The Heart of PM

  Communicating operators lie at the heart of the PM parallelisation model. They are designed to provide a compromise between the direct access to global data structures (particularly arrays) offered by approaches such as Partitioned Global Address Space  and the straightforward synchronisation provided by Communicating Sequential Processes . I n common with most data-parallel languages, PM contains a parallel version of the for statement which runs all invocations of its enclosed statement list concurrently: for element1 in array1, element2 in array2 do element2=process_element(element1) endfor Most real models will require some interaction between adjacent array elements. In PM this is achieved by using either a local or a global communicating operator. The global operator @v returns an array whose elements comprise the values of loop-local variable v in each invocation of the enclosing for statement. The neighbourhood operator v@{nbd ...
Post a Comment
Read more