
\documentclass{article}

\begin{document}

\title{2-D Potts Model Simulation using Metropolis Heatbath
in HPJava}

\author{Yuhong Wen, Bryan Carpenter}

\maketitle

\section{Introduction}

The initial 2-D Potts Model Simulation using Metropolis
Heatbath has been implemented in Fortran 90 and HPF in NPAC. We will
introduce its implementatin in HPJava here.

The code sets up a 2-D Q-State Potts model, and evolves it using the
classic Metropolis heatbath algorithm. In this model each site may take on
one of Q states, and each site interacts with its neighbour if both
states are equal.

\section{Implementation}

In terms of communications in the parallelised HPF code, the main cost comes
from the use of the CSHIFT() intrinsic functions in directions where the
data has been distributed, since the 2-D array of sites is block distributed
in one direction. Data transfer is near-neighbour between sites.

In HPJava implementation, we use the new feature of communication
{\tt writeHalo},
which only copies the necessary neighbour elements between the processors,
instead of moving the whole data along the sepecific dimention of the array.
This will dramaticly reduce the communication cost and memory copy operations
when updating the energy states.

Another difference in the HPJava implementation is that we can provide an
efficient and parallelized implementation of the random number generator
locally in each processor in HPJava. The locally controlled operations and
local temporary values will  particularly be of benefit to spin physics
codes such as this Potts model.

{\em [This example does not respect the standard variable usage rules
in the case of random number state {\tt xsubi}.  {\tt xsubi} is not a
distributed array, but it has different values on different processes.
Whether this should be allowed is an open question.  Incidentally,
we do not endorse the specific random number generator or mode
of initialization for accurate simulations!]}

\newpage
\scriptsize
\begin{verbatim}
class Potts {
  int wlo [[]] = {1, 1}, whi [[]]]] = {1, 1} ;        // ghost width for 'writeHalo'
  int mode [[]] = {CYCL, CYCL} ;          // boundary conds for 'writeHalo'

  void main(int argc, char **argv) {

    int nprocs = AdlibNProcs() ;

    Procs1 p = new Procs1(nprocs) ;

    on(p) {
        const int w [2] = {1, 1} ;      // ghost area for 'writeHalo'
        int mode[] = {CYCL, CYCL} ;     // boundary conds for 'writeHalo'

        int m = 64;
        int n = 64;

        Range x = new Range(m) ;
        Range y = new BlockRange(n, p.dim(0)) ;

        int q = 4;

        float cup = -4.0;

        int latt[[,]] = new int [[x, y]] ;
        int latt_old[[,]] = new int [[x, y]] ;

        char result1[[,]] = new char [[x, y]] ;
        char result2[[,]] = new char [[x, y]] ;
        char result3[[,]] = new char [[x, y]] ;

        int itmp[[,]] = new int [[x, y]] ;


        int nmeasure = 20;
        int niter = 10;

        int xsubi [3] = new int [3] ;

        float t1, t2;

//  Set up a random number generator.
        setup_random( xsubi, 12345 );

//  Initialize the latt array.
        over(x) {
            over(y) {
                latt[[x, y]] = int( erand48(xsubi) * q ) ;
            }
        }

        for (int imeasure = 1; imeasure<=nmeasure; imeasure++) {

//  Use ghost area to get all four neighbor values.
            Adlib.copy(latt_old, latt) ;
            Adlib.writeHalo(latt_old, w, w, mode) ;

            over (x) {
                over (y) {
                    int latt_val = latt[[x, y]] ;
                    itmp[[x, y]] = latt_val == latt_old[[x >> 1, y]] +
                                   latt_val == latt_old[[x << 1, y]] +
                                   latt_val == latt_old[[x, y >> 1]] +
                                   latt_val == latt_old[[x, y << 1]] ;
                }
            }

            Adlib.gprintf("\t# Bonds = %d\n",sum( itmp ));

            t1 = (float)clock()/CLOCKS_PER_SEC;
            for (int iter=1; iter<=niter; iter++) {

              for (int parity = 0; parity < 2; parity ++) {

                    Adlib.copy(latt_old, latt) ;
                    Adlib.writeHalo(latt_old, w, w, mode) ;

                    over (x) {
                      over (y | (parity + x) % 2 : : 2) {
                        int latt_val = latt[[x, y]] ;
                        int itmp_val = (latt[[x, y]] +
                                       int(erand48(xsubi) * q - 1) + 1) % q ;
                        float  p  = (latt_val == latt_old[[x >> 1, y]] +
                                     latt_val == latt_old[[x << 1, y]] +
                                     latt_val == latt_old[[x, y >> 1]] +
                                     latt_val == latt_old[[x, y << 1]] -
                                     itmp_val == latt_old[[x >> 1, y]] -
                                     itmp_val == latt_old[[x << 1, y]] -
                                     itmp_val == latt_old[[x, y >> 1]] -
                                     itmp_val == latt_old[[x, y << 1]]) * cup ;
                        if (erand48(xsubi) < exp( p ) )
                                latt[[x, y]] = itmp_val ;
                      }
                    }

               }

            }
            t2 = (float)clock()/CLOCKS_PER_SEC;
            Adlib.gprintf("Time: %f\n", t2 - t1);
        }

//  Output final value of `latt'
        FILE* ifile = Adlib.gfopen("out.ppm", "w") ;

        Adlib.gfprintf(ifile, "P6\n");
        Adlib.gfprintf(ifile, "# 2d Potts Model\n");
        Adlib.gfprintf(ifile, "%12d%12d\n",m,n);
        Adlib.gfprintf(ifile, "%8d\n",127);

        over(x) {
            over(y) {
                switch (latt[[x, y]]) {
                  case 0 :
                        result1[[x, y]] = 127;
                        result2[[x, y]] = 0;
                        result3[[x, y]] = 0;
                        break;
                  case 1 :
                        result1[[x, y]] = 0;
                        result2[[x, y]] = 127;
                        result3[[x, y]] = 0;
                        break;
                  case 2 :
                        result1[[x, y]] = 0;
                        result2[[x, y]] = 0;
                        result3[[x, y]] = 127;
                        break;
                  case 3 :
                        result1[[x, y]] = 127;
                        result2[[x, y]] = 127;
                        result3[[x, y]] = 0;
                        break;
                  case 4 :
                        result1[[x, y]] = 127;
                        result2[[x, y]] = 0;
                        result3[[x, y]] = 127;
                        break;
                  case 5 :
                        result1[[x, y]] = 0;
                        result2[[x, y]] = 127;
                        result3[[x, y]] = 127;
                        break;
                  case 6 :
                        result1[[x, y]] = 127;
                        result2[[x, y]] = 127;
                        result3[[x, y]] = 127;
                        break;
                  default :
                        result1[[x, y]] = 0;
                        result2[[x, y]] = 0;
                        result3[[x, y]] = 0;
                        break;
                }
            }
        }

        Adlib.afprintf(ifile, "%c%c%c",result1,result2,result3);

        Adlib.gfclose(ifile);

    }

  }
}


void setup_random(int xsubi [] , int iseed) {

  over(id) {
    int procs = id ;
    for (int i=0; i<3; i++) {
        xsubi[i] = (procs * 4096 + 889 + iseed) % 24025;
        procs *= 3;
    }
  }
}

\end{verbatim}
\normalsize

\newpage
\section{Performance}

The most time consuming part in the kernel loop when updating the energy
state of each site has been reduced by using the ghost area
communication technology {\tt writeHalo} in the HPJava implementation.
This will result in much better performance compared with the HPF and
F90 version.

(To be tested)

\end{document}

