Adding description of 2d wave equation discretization, and simplifying example code.
authorIsaac Dooley <idooley2@illinois.edu>
Thu, 16 Oct 2008 21:02:23 +0000 (21:02 +0000)
committerIsaac Dooley <idooley2@illinois.edu>
Thu, 16 Oct 2008 21:02:23 +0000 (21:02 +0000)
examples/charm++/wave2d/paper.pdf [new file with mode: 0644]
examples/charm++/wave2d/paper.tex [new file with mode: 0644]
examples/charm++/wave2d/runserver.sh
examples/charm++/wave2d/wave2d.C

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, 
+\begin{equation}
+\frac{\partial^2p}{\partial p^2}=c^2 \nabla^2 p ,\label{PDE}
+\end{equation}
+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:
+
+\begin{equation}
+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}
+\end{equation}
+
+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, 
+
+\begin{equation}
+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} .
+\end{equation}
+
+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}
+
index c49cd0f7c249c1fdd11cc077ac1d489d1cdbcde0..866c684c66cd466042cbde5401aadad450a82cde 100755 (executable)
@@ -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)
 /*readonly*/ CProxy_Main mainProxy;
 /*readonly*/ CProxy_Wave arrayProxy;
 
-#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;
+      int xcenter = radius + rand() % (TotalDataWidth - 2*radius);
+      int ycenter = radius + rand() % (TotalDataHeight - 2*radius);
+      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;
   }