Removing memory leak in kNeighbor benchmark
[charm.git] / examples / charm++ / wave2d / wave2d.C
1 #include "liveViz.h"
2 #include "wave2d.decl.h"
3
4 // This program solves the 2-d wave equation over a grid, displaying pretty results through liveViz
5 // The discretization used below is described in the accompanying paper.pdf
6 // Author: Isaac Dooley 2008
7
8 /*readonly*/ CProxy_Main mainProxy;
9 /*readonly*/ CProxy_Wave arrayProxy;
10
11 #define TotalDataWidth  800
12 #define TotalDataHeight 700
13 #define chareArrayWidth  4
14 #define chareArrayHeight  3
15 #define total_iterations 5000
16 #define numInitialPertubations 5
17
18 #define mod(a,b)  (((a)+b)%b)
19
20 enum { left=0, right, up, down };
21
22 class Main : public CBase_Main
23 {
24 public:
25   int iteration;
26   int count;
27
28   Main(CkArgMsg* m) {
29     iteration = 0;
30     count = 0;
31     mainProxy = thisProxy; // store the main proxy
32     
33     CkPrintf("Running wave2d on %d processors\n", CkNumPes());
34
35     // Create new array of worker chares
36     arrayProxy = CProxy_Wave::ckNew(chareArrayWidth, chareArrayHeight);
37
38     // setup liveviz
39     CkCallback c(CkIndex_Wave::requestNextFrame(0),arrayProxy);
40     liveVizConfig cfg(liveVizConfig::pix_color,true);
41     liveVizInit(cfg,arrayProxy,c);
42
43     //Start the computation
44     arrayProxy.begin_iteration();
45   }
46
47   // Each worker calls this method
48   void iterationCompleted() {
49     count++;
50     if(count == chareArrayWidth*chareArrayHeight){
51       if (iteration == total_iterations) {
52         CkPrintf("Program Done!\n");
53         CkExit();
54       } else { 
55         // Start the next iteration
56         count = 0;
57         iteration++;
58         if(iteration % 20 == 0) CkPrintf("Completed %d iterations\n", iteration);    
59         arrayProxy.begin_iteration();
60       }
61     }
62   }
63   
64   
65 };
66
67 class Wave: public CBase_Wave {
68 public:
69   int messages_due;
70   int mywidth;
71   int myheight;
72
73   double *pressure_old;  // time t-1
74   double *pressure; // time t
75   double *pressure_new;  // time t+1
76
77   double *buffers[4];
78
79   // Constructor, initialize values
80   Wave() {
81
82     mywidth=TotalDataWidth / chareArrayWidth;
83     myheight= TotalDataHeight / chareArrayHeight;
84
85     pressure_new  = new double[mywidth*myheight];
86     pressure = new double[mywidth*myheight];
87     pressure_old  = new double[mywidth*myheight];
88
89     buffers[left] = new double[myheight];
90     buffers[right]= new double[myheight];
91     buffers[up]   = new double[mywidth];
92     buffers[down] = new double[mywidth];
93
94     messages_due = 4;
95
96     InitialConditions();
97   }
98
99
100   // Setup some Initial pressure pertubations for timesteps t-1 and t
101   void InitialConditions(){
102     srand(0); // Force the same random numbers to be used for each chare array element
103
104     for(int i=0;i<myheight*mywidth;i++)
105       pressure[i] = pressure_old[i] = 0.0;
106     
107     for(int s=0; s<numInitialPertubations; s++){    
108       // Determine where to place a circle within the interior of the 2-d domain
109       int radius = 20+rand() % 30;
110       int xcenter = radius + rand() % (TotalDataWidth - 2*radius);
111       int ycenter = radius + rand() % (TotalDataHeight - 2*radius);
112       // Draw the circle
113       for(int i=0;i<myheight;i++){
114         for(int j=0; j<mywidth; j++){
115           int globalx = thisIndex.x*mywidth + j; // The coordinate in the global data array (not just in this chare's portion)
116           int globaly = thisIndex.y*myheight + i;
117           double distanceToCenter = sqrt((globalx-xcenter)*(globalx-xcenter) + (globaly-ycenter)*(globaly-ycenter));
118           if (distanceToCenter < radius) {
119             double rscaled = (distanceToCenter/radius)*3.0*3.14159/2.0; // ranges from 0 to 3pi/2 
120             double t = 700.0 * cos(rscaled) ; // Range won't exceed -700 to 700
121             pressure[i*mywidth+j] = pressure_old[i*mywidth+j] = t;
122           }
123         }                                               
124       }
125     }
126   }
127
128   Wave(CkMigrateMessage* m) { }
129
130   ~Wave() { }
131
132   void begin_iteration(void) {
133
134     double *top_edge = &pressure[0];
135     double *bottom_edge = &pressure[(myheight-1)*mywidth];
136
137     double *left_edge = new double[myheight];
138     double *right_edge = new double[myheight];
139     for(int i=0;i<myheight;++i){
140       left_edge[i] = pressure[i*mywidth];
141       right_edge[i] = pressure[i*mywidth + mywidth-1];
142     }
143
144     // Send my left edge
145     thisProxy(mod(thisIndex.x-1, chareArrayWidth), thisIndex.y).recvGhosts(right, myheight, left_edge);
146     // Send my right edge
147     thisProxy(mod(thisIndex.x+1, chareArrayWidth), thisIndex.y).recvGhosts(left, myheight, right_edge);
148     // Send my top edge
149     thisProxy(thisIndex.x, mod(thisIndex.y-1, chareArrayHeight)).recvGhosts(down, mywidth, top_edge);
150     // Send my bottom edge
151     thisProxy(thisIndex.x, mod(thisIndex.y+1, chareArrayHeight)).recvGhosts(up, mywidth, bottom_edge);
152
153     delete [] right_edge;
154     delete [] left_edge;
155   }
156   
157   void recvGhosts(int whichSide, int size, double ghost_values[]) {
158     for(int i=0;i<size;++i)
159       buffers[whichSide][i] = ghost_values[i];   
160     check_and_compute();
161   }
162
163   void check_and_compute() {
164     if (--messages_due == 0) {
165
166       // Compute the new values based on the current and previous step values
167
168       for(int i=0;i<myheight;++i){
169         for(int j=0;j<mywidth;++j){
170
171           // Current time's pressures for neighboring array locations
172           double L = (j==0          ? buffers[left][i]  : pressure[i*mywidth+j-1] );
173           double R = (j==mywidth-1  ? buffers[right][i] : pressure[i*mywidth+j+1] );
174           double U = (i==0          ? buffers[up][j]    : pressure[(i-1)*mywidth+j] );
175           double D = (i==myheight-1 ? buffers[down][j]  : pressure[(i+1)*mywidth+j] );
176
177           // Current time's pressure for this array location
178           double curr = pressure[i*mywidth+j];
179
180           // Previous time's pressure for this array location
181           double old  = pressure_old[i*mywidth+j];
182
183           // Compute the future time's pressure for this array location
184           pressure_new[i*mywidth+j] = 0.4*0.4*(L+R+U+D - 4.0*curr)-old+2.0*curr;
185
186         }
187       }
188                 
189       // Advance to next step by shifting the data back one step in time
190       double *tmp = pressure_old;
191       pressure_old = pressure;
192       pressure = pressure_new;
193       pressure_new = tmp;
194
195       messages_due = 4;
196       mainProxy.iterationCompleted();
197     }
198   }
199
200
201   // provide my portion of the image to the graphical liveViz client                           
202   void requestNextFrame(liveVizRequestMsg *m){
203
204     // Draw my part of the image, plus a nice 1px border along my right/bottom boundary
205     int sx=thisIndex.x*mywidth; // where my portion of the image is located
206     int sy=thisIndex.y*myheight;
207     int w=mywidth; // Size of my rectangular portion of the image
208     int h=myheight;
209     
210     // set the output pixel values for my rectangle
211     // Each RGB component is a char which can have 256 possible values.
212     unsigned char *intensity= new unsigned char[3*w*h];
213     for(int i=0;i<myheight;++i){
214       for(int j=0;j<mywidth;++j){
215
216         double p = pressure[i*mywidth+j];
217         if(p > 255.0) p = 255.0;    // Keep values in valid range
218         if(p < -255.0) p = -255.0;  // Keep values in valid range
219                 
220         if(p > 0) { // Positive values are red
221           intensity[3*(i*w+j)+0] = 255; // RED component
222           intensity[3*(i*w+j)+1] = 255-p; // GREEN component
223           intensity[3*(i*w+j)+2] = 255-p; // BLUE component
224         } else { // Negative values are blue
225           intensity[3*(i*w+j)+0] = 255+p; // RED component
226           intensity[3*(i*w+j)+1] = 255+p; // GREEN component
227           intensity[3*(i*w+j)+2] = 255; // BLUE component
228         }
229         
230       }
231     }
232     
233     // Draw a green border on right and bottom of this chare array's pixel buffer. 
234     // This will overwrite some pressure values at these pixels.
235     for(int i=0;i<h;++i){
236       intensity[3*(i*w+w-1)+0] = 0;     // RED component
237       intensity[3*(i*w+w-1)+1] = 255;   // GREEN component
238       intensity[3*(i*w+w-1)+2] = 0;     // BLUE component
239     }
240     for(int i=0;i<w;++i){
241       intensity[3*((h-1)*w+i)+0] = 0;   // RED component
242       intensity[3*((h-1)*w+i)+1] = 255; // GREEN component
243       intensity[3*((h-1)*w+i)+2] = 0;   // BLUE component
244     }
245
246     liveVizDeposit(m, sx,sy, w,h, intensity, this);
247     delete[] intensity;
248
249   }
250   
251 };
252
253 #include "wave2d.def.h"