redblack2d.c

LPARX is a C++ class library for for implementing portable scientific applications on distributed me

//
// Perform Jacobi iterations for a two dimensional Laplace's equation
// This is intended to be an example of how to use LPARX.  If you ever use
// Jacobi iterations to solve an elliptic problem, you will be severly punished
//

#include "redblack2d.h"

/*
  Each LPARX grid has two attributes:

	1. the codomain, the base type. e.g. double
	2. the rank, the number of dimensions

  LPARX codomains can be of any type.
  The current implementation supports ranks 1 through 4.

  The name of a Grid type is formed as follows

	GridX(T)

  where X is the rank (an int in the range 1..4) and T is a type.

  Macros for types Grid1, Grid2, Grid3, and Grid4 are defined in
  the standard header file "lparx.h", which must be included
  in every LPARX module (this include is found in "hdr.h".)
  We use Grid macros to declare and implement specific instances of Grid
*/

implement_Grid2(double)

/*
  For each grid type we must have a declaration
  and corresponding implementation.
  Here we implement a 2-D Grid class of double
  The corresponding declaration appears before it in the header file "hdr.h"
  Each module using a particular Grid type must declare it, i.e. in "hdr.h"
  The implementation must appear in exactly one .C module (generally main)
  This is a standard C++ way of emulating templates. 
  Implementations look like function calls, but in fact they invoke macros
  Because type double is a built-in type, we do no have to 
  define any pack-unpack routines for double.
  However, if the type were a user-defined class or structure,
  then we would need to provide hints to LPARX about how to
  gather and scatter the data type. This is relatively easy
  to do and is described in the User's Guide.
*/

//
// Forward and external procedure declarations
//

void initData(XArray1(Grid2(double))& U);
double relax(XArray1(Grid2(double))& U, const double& h);

#define OUTPUT(X) { if (mpMyID() == 0) cout << X << flush; }

//
// Beginning of main() -- define defaults for the problem
//

#define M 40
#define N 50
#define MAX_ITERATIONS 10
#define DELTA_RESIDUAL 1.0e-2

int main(const int argc, char **argv)
  {
   /*
     Define the problem domain to be [1:M, 1:N]
     A Region2 is a first-class object that describes an index domain.
   */

   Region2 Interior(1, 1, M, N);
   const double h = 1.0;

   /*
     Call mpNodes to determine the number of processors. 
   */

   int P = mpNodes();

   /*
     We will partition the problem so that each processor gets only
     one region. Function UniformPartition() divides the domain into
     P regions in a block decomposition, each region approximately
     the same size. UniformPartition() returns a pointer to a list of
     these regions. For more information on UniformPartition, consult
     the LPARX documentation.
   */

   Region2 *Partition = UniformPartition(Interior, P);

   /*
     We need a partition with room to hold one
     ghost cell value in each direction.
     Function grow() creates this new list of regions.
   */

   Region2 *Ghosts = grow(Partition, P, 1);

   /*
     Declare the XArray of grids.  The XArray1 constructor declares
     a one-dimensional array of placeholders for grids. The XAlloc()
     procedure instantiates the Grid storage for each grid in an XArray. 
   */

   XArray1(Grid2(double)) U;
   XAlloc(U, P, Ghosts);

   /* Set up the initial conditions */

   initData(U);

   /*
     Perform the iterations
   */

   for (int i = 0; i < MAX_ITERATIONS; i++)
     {
      const double dU = relax(U, h);
      OUTPUT("Delta at iteration " << i << ": " << dU << endl);
      if (dU < DELTA_RESIDUAL) break;
     }

   delete [] Partition;
   delete [] Ghosts;
  }