2929#include <time.h>
3030
3131#include <paraconf.h>
32+ #include <pdi.h>
3233
33- /// size of the local data as [HEIGHT, WIDTH] including ghosts & boundary constants
34+ // size of the local data as [HEIGHT, WIDTH] including the number of ghost layers
35+ // for communications or boundary conditions
3436int dsize [2 ];
3537
36- /// 2D size of the process grid as [HEIGHT, WIDTH]
38+ // 2D size of the process grid as [HEIGHT, WIDTH]
3739int psize [2 ];
3840
39- /// 2D rank of the local process in the process grid as [YY, XX]
41+ // 2D rank of the local process in the process grid as [YY, XX]
4042int pcoord [2 ];
4143
42- /// the alpha coefficient used in the computation
44+ // the alpha coefficient used in the computation
4345double alpha ;
4446
45- /** Initialize the data all to 0 except for the left border (XX==0) initialized to 1 million
47+ double L = 1.0 ;
48+ // definition of the source
49+ // the source corresponds to a disk of an uniform value
50+ // source1: center=(0.4,0.4), radius=0.2 and value=100
51+ double source1 [4 ]= {0.4 , 0.4 , 0.2 , 100 };
52+ // source2: center=(0.8,0.7), radius=0.1 and value=200
53+ double source2 [4 ]= {0.7 , 0.8 , 0.1 , 200 };
54+ // the order of the coordinates of the center (XX,YY) is inverted in the vector
55+
56+ /** Initialize all the data to 0, with the exception of a given cell
57+ * whose center (cpos_x,cpos_y) is inside of the disks
58+ * defined by source1 or source2
4659 * \param[out] dat the local data to initialize
4760 */
4861void init (double dat [dsize [0 ]][dsize [1 ]])
4962{
5063 for (int yy = 0 ; yy < dsize [0 ]; ++ yy ) for (int xx = 0 ; xx < dsize [1 ]; ++ xx ) dat [yy ][xx ] = 0 ;
51- if ( pcoord [1 ] == 0 ) for (int yy = 0 ; yy < dsize [0 ]; ++ yy ) dat [yy ][0 ] = 1000000 ;
64+ double dy = L / ((dsize [0 ]- 2 ) * psize [0 ]) ;
65+ double dx = L / ((dsize [1 ]- 2 ) * psize [1 ]) ;
66+
67+ double cpos_x ,cpos_y ;
68+ double square_dist1 , square_dist2 ;
69+ for (int yy = 0 ; yy < dsize [0 ];++ yy ) {
70+ cpos_y = (yy + pcoord [0 ]* (dsize [0 ]- 2 ))* dy - 0.5 * dy ;
71+ for (int xx = 0 ; xx < dsize [1 ];++ xx ) {
72+ cpos_x = (xx + pcoord [1 ]* (dsize [1 ]- 2 ))* dx - 0.5 * dx ;
73+ square_dist1 = ( cpos_y - source1 [0 ] ) * ( cpos_y - source1 [0 ] )
74+ + ( cpos_x - source1 [1 ] ) * ( cpos_x - source1 [1 ] );
75+ if (square_dist1 <= source1 [2 ] * source1 [2 ]) {
76+ dat [yy ][xx ] = source1 [3 ];
77+ }
78+ square_dist2 = ( cpos_y - source2 [0 ] ) * ( cpos_y - source2 [0 ] )
79+ + ( cpos_x - source2 [1 ] ) * ( cpos_x - source2 [1 ] );
80+ if (square_dist2 <= source2 [2 ] * source2 [2 ]) {
81+ dat [yy ][xx ] = source2 [3 ];
82+ }
83+ }
84+ }
5285}
5386
5487/** Compute the values at the next time-step based on the values at the current time-step
@@ -58,24 +91,18 @@ void init(double dat[dsize[0]][dsize[1]])
5891void iter (double cur [dsize [0 ]][dsize [1 ]], double next [dsize [0 ]][dsize [1 ]])
5992{
6093 int xx , yy ;
61- for (xx = 0 ; xx < dsize [1 ]; ++ xx ) next [0 ][xx ] = cur [0 ][xx ];
6294 for (yy = 1 ; yy < dsize [0 ]- 1 ; ++ yy ) {
63- next [yy ][0 ] = cur [yy ][0 ];
6495 for (xx = 1 ; xx < dsize [1 ]- 1 ; ++ xx ) {
65- next [yy ][xx ] =
66- (1. - 4. * alpha ) * cur [yy ][xx ]
67- + alpha * ( cur [yy ][xx - 1 ]
68- + cur [yy ][xx + 1 ]
69- + cur [yy - 1 ][xx ]
70- + cur [yy + 1 ][xx ]
71- );
96+ next [yy ][xx ] = (1. - 4. * alpha ) * cur [yy ][xx ]
97+ + alpha * ( cur [yy ][xx - 1 ]
98+ + cur [yy ][xx + 1 ]
99+ + cur [yy - 1 ][xx ]
100+ + cur [yy + 1 ][xx ]);
72101 }
73- next [yy ][dsize [1 ]- 1 ] = cur [yy ][dsize [1 ]- 1 ];
74102 }
75- for (xx = 0 ; xx < dsize [1 ]; ++ xx ) next [dsize [0 ]- 1 ][xx ] = cur [dsize [0 ]- 1 ][xx ];
76103}
77104
78- /** Exchanges ghost values with neighbours
105+ /** Exchange ghost values with neighbours
79106 * \param[in] cart_comm the MPI communicator with all processes organized in a 2D Cartesian grid
80107 * \param[in] cur the local data at the current time-step whose ghosts need exchanging
81108 */
@@ -102,8 +129,8 @@ void exchange(MPI_Comm cart_comm, double cur[dsize[0]][dsize[1]])
102129
103130 // send up
104131 MPI_Cart_shift (cart_comm , 0 , -1 , & rank_source , & rank_dest );
105- MPI_Sendrecv (& cur [1 ][1 ], 1 , row , rank_dest , 100 , // send column after ghost
106- & cur [dsize [0 ]- 1 ][1 ], 1 , row , rank_source , 100 , // receive last column (ghost)
132+ MPI_Sendrecv (& cur [1 ][1 ], 1 , row , rank_dest , 100 , // send row after ghost
133+ & cur [dsize [0 ]- 1 ][1 ], 1 , row , rank_source , 100 , // receive last row (ghost)
107134 cart_comm , & status );
108135
109136 // send to the right
@@ -114,7 +141,7 @@ void exchange(MPI_Comm cart_comm, double cur[dsize[0]][dsize[1]])
114141
115142 // send to the left
116143 MPI_Cart_shift (cart_comm , 1 , -1 , & rank_source , & rank_dest );
117- MPI_Sendrecv (& cur [1 ][1 ], 1 , column , rank_dest , 100 , // send column after ghost
144+ MPI_Sendrecv (& cur [1 ][1 ], 1 , column , rank_dest , 100 , // send column after ghost
118145 & cur [1 ][dsize [1 ]- 1 ], 1 , column , rank_source , 100 , // receive last column (ghost)
119146 cart_comm , & status );
120147}
@@ -138,8 +165,9 @@ int main( int argc, char* argv[] )
138165 // load the alpha parameter
139166 PC_double (PC_get (conf , ".alpha" ), & alpha );
140167
141- // load the global data-size
142168 int global_size [2 ];
169+ // load the global data-size
170+ // you can use paraconf to read some parameters from the yml config file
143171 PC_int (PC_get (conf , ".global_size.height" ), & longval ); global_size [0 ] = longval ;
144172 PC_int (PC_get (conf , ".global_size.width" ), & longval ); global_size [1 ] = longval ;
145173
@@ -152,12 +180,12 @@ int main( int argc, char* argv[] )
152180 assert (global_size [1 ]%psize [1 ]== 0 );
153181 assert (psize [1 ]* psize [0 ] == psize_1d );
154182
155- // compute the local data-size with space for ghosts and boundary constants
183+ // compute the local data-size (the number of ghost layers is 2 for each coordinate)
156184 dsize [0 ] = global_size [0 ]/psize [0 ] + 2 ;
157185 dsize [1 ] = global_size [1 ]/psize [1 ] + 2 ;
158186
159187 // create a 2D Cartesian MPI communicator & get our coordinate (rank) in it
160- int cart_period [2 ] = { 0 , 0 };
188+ int cart_period [2 ] = { 1 , 1 };
161189 MPI_Comm cart_comm ; MPI_Cart_create (main_comm , 2 , psize , cart_period , 1 , & cart_comm );
162190 MPI_Cart_coords (cart_comm , pcoord_1d , 2 , pcoord );
163191
@@ -173,6 +201,7 @@ int main( int argc, char* argv[] )
173201
174202 // the main loop
175203 for (; ii < 10 ; ++ ii ) {
204+
176205 // compute the values for the next iteration
177206 iter (cur , next );
178207
0 commit comments