author Isaac Dooley Thu, 16 Oct 2008 21:02:23 +0000 (21:02 +0000) committer Isaac Dooley Thu, 16 Oct 2008 21:02:23 +0000 (21:02 +0000)
 examples/charm++/wave2d/paper.pdf [new file with mode: 0644] patch | blob examples/charm++/wave2d/paper.tex [new file with mode: 0644] patch | blob examples/charm++/wave2d/runserver.sh patch | blob | history examples/charm++/wave2d/wave2d.C patch | blob | history

diff --git a/examples/charm++/wave2d/paper.pdf b/examples/charm++/wave2d/paper.pdf
new file mode 100644 (file)
index 0000000..f62ca29
Binary files /dev/null and b/examples/charm++/wave2d/paper.pdf differ
diff --git a/examples/charm++/wave2d/paper.tex b/examples/charm++/wave2d/paper.tex
new file mode 100644 (file)
index 0000000..38bfea0
--- /dev/null
@@ -0,0 +1,36 @@
+% Authors: Isaac Dooley, Sandhya Mangala, Laxmikant Kale, and Philippe Geubelle
+% Journal: Submitting to CSD
+
+
+\documentclass{article}
+
+%-------------------------------------------------------------------------
+\begin{document}
+
+\title{2-D Wave Equation}
+\author{ Isaac Dooley }
+
+
+
+\maketitle
+
+The wave equation can be stated in the following form,
+
+\frac{\partial^2p}{\partial p^2}=c^2 \nabla^2 p ,\label{PDE}
+
+where p is a physical measure such as a pressure, or the height of the water in a tub. This equation can be discretized over a 2d grid using a finite differencing scheme. Assume the pressures are located across a rectangular 2-D grid in x and y dimensions. We call the value of the pressure $p_{x,y}^{t}$ for time $t$ at location $\left(x,y\right)$ in the grid. One solution to the differential equation (\ref{PDE}) is the following:
+
+
+p_{x,y}^{t+1} + p_{x,y}^{t-1} -2 p_{x,y}^{t} = c^2 \left(p_{x+1,y}^{t}+p_{x-1,y}^{t}+p_{x,y+1}^{t}+p_{x,y-1}^{t}-4p_{x,y}^{t}\right).\label{DISC}
+
+
+The left hand side term $\frac{\partial^2p}{\partial p^2}$ is represented in terms of the values of $p$ at three consecutive timesteps at grid location $\left(x,y\right)$. The right hand side terms are discretized using the $p$ values for the 4 locations adjacent to $\left(x,y\right)$ in the grid. One can simply solve for $p_{x,y}^{t+1}$ in equation (\ref{DISC}) to produce an update rule,
+
+
+p_{x,y}^{t+1} = c^2 \left(p_{x+1,y}^{t}+p_{x-1,y}^{t}+p_{x,y+1}^{t}+p_{x,y-1}^{t}-4p_{x,y}^{t}\right)- p_{x,y}^{t-1} +2 p_{x,y}^{t} .
+
+
+This update rule specifies how the $p$ value for each location in the grid should be computed, using $p$ values from the two previous time steps. This update rule is easy to code in an application. The application simply maintains two timesteps' worth of $p$ grids, and using these, the next timestep's $p$ grid can be computed.
+
+\end{document}
+
@@ -1,5 +1,5 @@
#!/bin/sh
echo "Starting up server:"
-./charmrun +p2 ./wave2d ++server ++server-port 1234
+./charmrun +p4 ./wave2d ++server ++server-port 1234
echo "Server Exited"

index 2c712dc93c2971915e3bd9915f91482c3685e572..9c61cf19c98b9d78d36484f0acfdede9fb755d11 100644 (file)

-#define TotalDataWidth  600
-#define TotalDataHeight 600
+#define TotalDataWidth  800
+#define TotalDataHeight 800

#define chareArrayWidth  4
#define chareArrayHeight  4
+#define num_chares ((chareArrayWidth)*(chareArrayHeight))
+
+#define total_iterations 5000
+
+#define numInitialPertubations 5
+

// A modulo operator that works for a==-1
#define wrap(a,b)  (((a)+b)%b)
class Main : public CBase_Main
{
public:
-       int recieve_count;
-       int num_chares;
-       int iterations;
-       int total_iterations;
-       double startTime;
-
-       Main(CkArgMsg* m) {
+  int recieve_count;
+  int iterations;

-               srand(0);
+  Main(CkArgMsg* m) {

-               // set iteration counter to zero
-               iterations=0;
+    // set iteration counter to zero
+    iterations=0;

-               // store the main proxy
-               mainProxy = thisProxy;
+    // store the main proxy
+    mainProxy = thisProxy;

-               // print info
-               CkPrintf("Running wave2d on %d processors\n", CkNumPes());
+    // print info
+    CkPrintf("Running wave2d on %d processors\n", CkNumPes());

-               total_iterations = 5000;
+    // Create new array of worker chares
+    arrayProxy = CProxy_Wave::ckNew(chareArrayWidth, chareArrayHeight);

-               // save the total number of worker chares we have in this simulation
-               num_chares = chareArrayHeight * chareArrayWidth;
-
-               // Create new array of worker chares
-               arrayProxy = CProxy_Wave::ckNew(chareArrayWidth, chareArrayHeight);
-
-               // setup liveviz
-               CkCallback c(CkIndex_Wave::requestNextFrame(0),arrayProxy);
-               liveVizConfig cfg(liveVizConfig::pix_color,true);
+    // setup liveviz
+    CkCallback c(CkIndex_Wave::requestNextFrame(0),arrayProxy);
+    liveVizConfig cfg(liveVizConfig::pix_color,true);

CkArrayID a; // unused???
-               liveVizInit(cfg,a,c);
-
-               //Start the computation
-               startTime = CmiWallTimer();
-               recieve_count = 0;
-               arrayProxy.begin_iteration();
-       }
-
+    liveVizInit(cfg,a,c);

+    //Start the computation
+    recieve_count = 0;
+    arrayProxy.begin_iteration();
+  }

-       // Each worker reports back to here when it completes an iteration
-       void report(int row, int col) {

-               recieve_count++;
-               double totaltime = CmiWallTimer() - startTime;

-               if (num_chares == recieve_count) {
+  // Each worker reports back to here when it completes an iteration
+  void report(int row, int col) {

-                       if (iterations == total_iterations) {
-                               CkPrintf("Completed %d iterations; last iteration time: %.6lf\n", iterations, totaltime);
-                               CkExit();
-                       } else {
+    recieve_count++;

-                         // Start the next iteration
-                         recieve_count=0;
-                         iterations++;
-                               CkPrintf("Completed %d iterations; iteration time: %.6lf\n", iterations, totaltime);
+    if (num_chares == recieve_count) {

-                         startTime = CmiWallTimer();
+      if (iterations == total_iterations) {
+       CkPrintf("Program Done!\n");
+       CkExit();
+      } else {

-                         for(int i=0;i<chareArrayWidth;i++){
-                           for(int j=0;j<chareArrayHeight;j++){
-                             arrayProxy(i,j).begin_iteration();
-                           }
-                         }
+       // Start the next iteration
+       recieve_count=0;
+       iterations++;
+       CkPrintf("Completed %d iterations\n", iterations);

-                       }
-               }
-       }
+       arrayProxy.begin_iteration();
+
+      }
+    }
+  }

};

class Wave: public CBase_Wave {
public:
-       int messages_due;
-       int mywidth;
-       int myheight;
+  int messages_due;
+  int mywidth;
+  int myheight;

-       double *temperature;
-       double *temperature_old;
+  double *pressure;
+  double *pressure_old;

-       double *buffer_left;
-       double *buffer_right;
-       double *buffer_up;
-       double *buffer_down;
+  double *buffer_left;
+  double *buffer_right;
+  double *buffer_up;
+  double *buffer_down;

-       // Constructor, initialize values
-       Wave() {
+  // Constructor, initialize values
+  Wave() {

-               mywidth=TotalDataWidth / chareArrayWidth;
-               myheight= TotalDataHeight / chareArrayHeight;
+    mywidth=TotalDataWidth / chareArrayWidth;
+    myheight= TotalDataHeight / chareArrayHeight;

-               temperature = new double[mywidth*myheight];
-               temperature_old = new double[mywidth*myheight];
+    pressure = new double[mywidth*myheight];
+    pressure_old = new double[mywidth*myheight];

-               buffer_left = new double[myheight];
+    buffer_left = new double[myheight];
buffer_right = new double[myheight];
buffer_up = new double[mywidth];
buffer_down = new double[mywidth];

-               messages_due = 4;
+    messages_due = 4;

-               InitialConditions();
-       }
-
-
-       // Setup some Initial conditions
-       void InitialConditions(){
+    InitialConditions();
+  }

-               int numInitialSpots = 5;
-               double spotScaleX = (double)TotalDataWidth / 800.0 ;
-               double spotScaleY = (double)TotalDataHeight / 800.0 ;
-               double spotScale = spotScaleX<spotScaleY ? spotScaleX : spotScaleY;

-               int xcenters[] = {50,300,420,580,380}; // X  coordinate of the initial pertubations
-               int ycenters[] = {50,300,420,580,511}; // Y coordinate of the initial pertubations
-               int rs[] = {10,30,10,20,15}; // radius of the initial pertubations
+  // Setup some Initial pressure pertubations
+  void InitialConditions(){
+    srand(0);

-    for(int i=0;i<myheight;i++){
-      for(int j=0; j<mywidth; j++){
-
-        temperature[i*mywidth+j] = 0.0;
-        temperature_old[i*mywidth+j] = 0.0;
-      }
+    for(int i=0;i<myheight*mywidth;i++){
+        pressure[i] = 0.0;
+        pressure_old[i] = 0.0;
}
-
-               for(int n=0;n<numInitialSpots; n++){
-
-                       double xcenter = (double)xcenters[n] * spotScale;
-                       double ycenter = (double)ycenters[n] * spotScale;
-                       double r = rs[n] * spotScale;
-
-                       for(int i=0;i<myheight;i++){
-                               for(int j=0; j<mywidth; j++){
-
-                                       int globalx = thisIndex.x*mywidth + j;
-                                       int globaly = thisIndex.y*myheight + i;
-
-                                       double thisr = sqrt((globalx-xcenter)*(globalx-xcenter) + (globaly-ycenter)*(globaly-ycenter));
-
-                                       if( thisr < r){
-                                               double rscaled = thisr / r; // ranges from 0 to 1
-                                               double rscaled2 = rscaled*3.0*3.14159/2.0; // ranges from 0 to 3pi/2
-
-                                               double t = 400.0 * cos(rscaled2) ; // Range between -1 to 1
-
-                                               temperature[i*mywidth+j] = t;
-                                               temperature_old[i*mywidth+j] = t;
-                                       }
-
-                               }
+
+    for(int s=0; s<numInitialPertubations; s++){
+      // Randomly place a circle within the 2-d data array (without wrapping around)
+      int radius = 20+rand() % 30;
+      for(int i=0;i<myheight;i++){
+       for(int j=0; j<mywidth; j++){
+         int globalx = thisIndex.x*mywidth + j; // The coordinate in the global data array (not just in this chare's portion)
+         int globaly = thisIndex.y*myheight + i;
+         double distanceToCenter = sqrt((globalx-xcenter)*(globalx-xcenter) + (globaly-ycenter)*(globaly-ycenter));
+         if (distanceToCenter < radius) {
+           double rscaled = distanceToCenter / radius; // ranges from 0 to 1
+           double rscaled2 = rscaled*3.0*3.14159/2.0; // ranges from 0 to 3pi/2
+           double t = 400.0 * cos(rscaled2) ; // Range not to exceed -400 to 400
+           pressure[i*mywidth+j] = t;
+           pressure_old[i*mywidth+j] = t;
+         }
+       }
}
-               }
+    }
+  }

-       }

+  Wave(CkMigrateMessage* m) {
+    CkAbort("Migration of this class is not supported yet. Write PUP function if migration is used\n");
+  }

-       Wave(CkMigrateMessage* m) {
-         CkAbort("Migration of this class is not supported yet. Write PUP function if migration is used\n");
-       }
+  ~Wave() {
+    delete [] pressure;
+    delete [] pressure_old;
+  }

-       ~Wave() {
-               delete [] temperature;
-               delete [] temperature_old;
-       }
+  void begin_iteration(void) {

-       void begin_iteration(void) {
-
-               double *left_edge = new double[myheight];
-               double *right_edge = new double[myheight];
-               double *top_edge = new double[mywidth];
-               double *bottom_edge = new double[mywidth];
-
-               for(int i=0;i<myheight;++i){
-                       left_edge[i] = temperature[i*mywidth];
-                       right_edge[i] = temperature[i*mywidth + mywidth-1];
-               }
-
-               for(int i=0;i<mywidth;++i){
-                       top_edge[i] = temperature[i];
-                       bottom_edge[i] = temperature[(myheight-1)*mywidth + i];
-               }
-
-               // Send my left edge
-               thisProxy(wrap(thisIndex.x-1,chareArrayWidth), thisIndex.y).ghostsFromRight(myheight, left_edge);
-               // Send my right edge
-               thisProxy(wrap(thisIndex.x+1,chareArrayWidth), thisIndex.y).ghostsFromLeft(myheight, right_edge);
-               // Send my top edge
-               thisProxy(thisIndex.x, wrap(thisIndex.y-1,chareArrayHeight)).ghostsFromBottom(mywidth, top_edge);
-               // Send my bottom edge
-               thisProxy(thisIndex.x, wrap(thisIndex.y+1,chareArrayHeight)).ghostsFromTop(mywidth, bottom_edge);
-
-               delete [] right_edge;
-               delete [] left_edge;
-               delete [] top_edge;
-               delete [] bottom_edge;
-       }
+    double *left_edge = new double[myheight];
+    double *right_edge = new double[myheight];
+    double *top_edge = new double[mywidth];
+    double *bottom_edge = new double[mywidth];

-       void ghostsFromRight(int width, double ghost_values[]) {
-               for(int i=0;i<width;++i){
-                       buffer_right[i] = ghost_values[i];
-               }
-               check_and_compute();
-       }
+    for(int i=0;i<myheight;++i){
+      left_edge[i] = pressure[i*mywidth];
+      right_edge[i] = pressure[i*mywidth + mywidth-1];
+    }

-       void ghostsFromLeft(int width, double ghost_values[]) {
-               for(int i=0;i<width;++i){
-                       buffer_left[i] = ghost_values[i];
-               }
-               check_and_compute();
-       }
+    for(int i=0;i<mywidth;++i){
+      top_edge[i] = pressure[i];
+      bottom_edge[i] = pressure[(myheight-1)*mywidth + i];
+    }

-       void ghostsFromBottom(int width, double ghost_values[]) {
-               for(int i=0;i<width;++i){
-                       buffer_down[i] = ghost_values[i];
-               }
-               check_and_compute();
-       }
+    // Send my left edge
+    thisProxy(wrap(thisIndex.x-1,chareArrayWidth), thisIndex.y).ghostsFromRight(myheight, left_edge);
+    // Send my right edge
+    thisProxy(wrap(thisIndex.x+1,chareArrayWidth), thisIndex.y).ghostsFromLeft(myheight, right_edge);
+    // Send my top edge
+    thisProxy(thisIndex.x, wrap(thisIndex.y-1,chareArrayHeight)).ghostsFromBottom(mywidth, top_edge);
+    // Send my bottom edge
+    thisProxy(thisIndex.x, wrap(thisIndex.y+1,chareArrayHeight)).ghostsFromTop(mywidth, bottom_edge);
+
+    delete [] right_edge;
+    delete [] left_edge;
+    delete [] top_edge;
+    delete [] bottom_edge;
+  }

-       void ghostsFromTop(int width, double ghost_values[]) {
-               for(int i=0;i<width;++i){
-                       buffer_up[i] = ghost_values[i];
-               }
-               check_and_compute();
-       }
+  void ghostsFromRight(int width, double ghost_values[]) {
+    for(int i=0;i<width;++i){
+      buffer_right[i] = ghost_values[i];
+    }
+    check_and_compute();
+  }

-       void check_and_compute() {
-               if (--messages_due == 0) {
-                       messages_due = 4;
-                       compute();
-                       mainProxy.report(thisIndex.x, thisIndex.y);
-               }
-       }
+  void ghostsFromLeft(int width, double ghost_values[]) {
+    for(int i=0;i<width;++i){
+      buffer_left[i] = ghost_values[i];
+    }
+    check_and_compute();
+  }

-       // We have recieved a message from all neighbors
-       void compute() {
+  void ghostsFromBottom(int width, double ghost_values[]) {
+    for(int i=0;i<width;++i){
+      buffer_down[i] = ghost_values[i];
+    }
+    check_and_compute();
+  }

-               // 2-d wave equation
-               // See http://www.mtnmath.com/whatrh/node66.html for description of this standard 2-d discretization of the wave equation
+  void ghostsFromTop(int width, double ghost_values[]) {
+    for(int i=0;i<width;++i){
+      buffer_up[i] = ghost_values[i];
+    }
+    check_and_compute();
+  }

+  void check_and_compute() {
+    if (--messages_due == 0) {

-    // Compute the new values based on the current and previous step values
+      // Compute the new values based on the current and previous step values

-               double temperature_new[mywidth*myheight];
-               for(int i=0;i<mywidth*myheight;i++)
-                       temperature_new[i] = 0.0;
+      double *pressure_new = new double[mywidth*myheight];

-               for(int i=0;i<myheight;++i){
-                       for(int j=0;j<mywidth;++j){
+      for(int i=0;i<myheight;++i){
+       for(int j=0;j<mywidth;++j){

-        // Current values for neighboring array elements
-                               double left  = (j==0          ? buffer_left[i]  : temperature[i*mywidth+j-1] );
-                               double right = (j==mywidth-1  ? buffer_right[i] : temperature[i*mywidth+j+1] );
-                               double up    = (i==0          ? buffer_up[j]    : temperature[(i-1)*mywidth+j] );
-                               double down  = (i==myheight-1 ? buffer_down[j]  : temperature[(i+1)*mywidth+j] );
+         // Current step's values for neighboring array elements
+         double left  = (j==0          ? buffer_left[i]  : pressure[i*mywidth+j-1] );
+         double right = (j==mywidth-1  ? buffer_right[i] : pressure[i*mywidth+j+1] );
+         double up    = (i==0          ? buffer_up[j]    : pressure[(i-1)*mywidth+j] );
+         double down  = (i==myheight-1 ? buffer_down[j]  : pressure[(i+1)*mywidth+j] );

-        // Current values for this array element
-                               double curr = temperature[i*mywidth+j];
+         // Current values for this array element
+         double curr = pressure[i*mywidth+j];

-        // Old value for this array element
-                               double old  = temperature_old[i*mywidth+j];
+         // Previous step's value for this array element
+         double old  = pressure_old[i*mywidth+j];

-        // The wave speed
-                               double c = 0.4;
+         // Wave speed
+         double c = 0.4;

-        // Compute the new value
-                               temperature_new[i*mywidth+j] = c*c*(left+right+up+down - 4.0*curr)-old+2.0*curr;
+         // Compute the new value
+         pressure_new[i*mywidth+j] = c*c*(left+right+up+down - 4.0*curr)-old+2.0*curr;

-        // Round any near-zero values to zero (avoid denorms)
-                               if(temperature_new[i*mywidth+j] < 0.0001 && temperature_new[i*mywidth+j] >  -0.0001)
-                                       temperature_new[i*mywidth+j] = 0.0;
+         // Round any near-zero values to zero (avoid denorms)
+         //      if(pressure_new[i*mywidth+j] < 0.0001 && pressure_new[i*mywidth+j] >  -0.0001)
+         //  pressure_new[i*mywidth+j] = 0.0;

-                       }
-               }
+       }
+      }

-               // Advance to next step by copying values to the arrays for the previous steps
-               for(int i=0;i<myheight;++i){
-                       for(int j=0;j<mywidth;++j){
-                               temperature_old[i*mywidth+j] = temperature[i*mywidth+j];
-                               temperature[i*mywidth+j] = temperature_new[i*mywidth+j];
-                       }
-               }
-
+      // Advance to next step by copying values to the arrays for the previous steps
+      for(int i=0;i<myheight;++i){
+       for(int j=0;j<mywidth;++j){
+         pressure_old[i*mywidth+j] = pressure[i*mywidth+j];
+         pressure[i*mywidth+j] = pressure_new[i*mywidth+j];
}
+      }
+
+      delete[] pressure_new;
+
+      messages_due = 4;
+      mainProxy.report(thisIndex.x, thisIndex.y);
+    }
+  }
+

// provide my portion of the image to the graphical liveViz client
@@ -332,18 +293,18 @@ public:
unsigned char *intensity= new unsigned char[3*w*h];
for(int i=0;i<myheight;++i){
for(int j=0;j<mywidth;++j){
-        double t = temperature[i*mywidth+j];
+        double t = pressure[i*mywidth+j];
if(t > 255.0){
t = 255.0;
} else if (t < -255.0){
t = -255.0;
}

-        if(t > 0) {
+        if(t > 0) { // Positive values are red
intensity[3*(i*w+j)+0] = 255; // RED component
intensity[3*(i*w+j)+1] = 255-t; // GREEN component
intensity[3*(i*w+j)+2] = 255-t; // BLUE component
-        } else {
+        } else { // Negative values are blue
intensity[3*(i*w+j)+0] = 255+t; // RED component
intensity[3*(i*w+j)+1] = 255+t; // GREEN component
intensity[3*(i*w+j)+2] = 255; // BLUE component
@@ -352,20 +313,19 @@ public:
}
}

-    // Draw border on right, overwrites data being plotted

+    // Draw a green border on right and bottom, overwrites data being plotted
for(int i=0;i<h;++i){
-      intensity[3*(i*w+w-1)+0] = 0; // RED component
+      intensity[3*(i*w+w-1)+0] = 0;     // RED component
intensity[3*(i*w+w-1)+1] = 255;   // GREEN component
-      intensity[3*(i*w+w-1)+2] = 0; // BLUE component
+      intensity[3*(i*w+w-1)+2] = 0;     // BLUE component
}
-
-    // Draw border on bottom, overwrites data being plotted
for(int i=0;i<w;++i){
-      intensity[3*((h-1)*w+i)+0] = 0; // RED component
-      intensity[3*((h-1)*w+i)+1] = 255;   // GREEN component
-      intensity[3*((h-1)*w+i)+2] = 0; // BLUE component
+      intensity[3*((h-1)*w+i)+0] = 0;   // RED component
+      intensity[3*((h-1)*w+i)+1] = 255; // GREEN component
+      intensity[3*((h-1)*w+i)+2] = 0;   // BLUE component
}
+
liveVizDeposit(m, sx,sy, w,h, intensity, this);
delete[] intensity;
}