/* tab:4 * * wator.sc - particle-like fish under current * * Author: Steve Lumetta * Version: 3 * Creation Date: Wed May 26 15:34:46 1993 * Filename: wator.sc * History: * SL 3 Wed Jun 2 11:19:37 1993 * Hacked force function into new version. * AK 2 Thu May 27 23:24:02 1993 * Completed version 2 based on steve's original version. * SL 1 Wed May 26 15:35:21 1993 * First written */ #include #include #include #include #include #define NFISH 100 /* number of fish */ #define DISPLAY_SIZE 256 /* pixel size of display window */ #define DISPLAY_RANGE 20 /* display range of fish space */ #define STEPS_PER_DISPLAY 10 /* time steps between display of fish */ #define T_FINAL 10.0 /* simulation end time */ #define GRAV_CONSTANT 1.0 /* proportionality constant of gravitational interaction */ #define MAX(X,Y) ((X) > (Y) ? (X) : (Y)) /* utility function */ /* This structure holds information for a single fish, including position, velocity, and mass. */ typedef struct { double x_pos, y_pos; double x_vel, y_vel; double mass; } fish_t; /* An atomic procedure which increments the number of fish present at a given gridpoint of the displayed ocean. */ void add_fish(int *ptr) { (*ptr)++; } /* Calculates a new ocean map from the fish positions, displaying the region [-DISPLAY_RANGE, DISPLAY_RANGE) x (same). */ static int *spread all_calculate_display(int num_fish, fish_t fish_list[]) { static int init = 0; static int *spread map; int i; fish_t *fish; int x_disp, y_disp; /* If this is the first time, allocate the map. */ if (!init) { init = 1; map = (int *spread) all_spread_malloc(DISPLAY_SIZE*DISPLAY_SIZE, sizeof(int)); } /* Set the new map to empty (colored blue), then synchronize before adding fish. */ for_my_1d(i, DISPLAY_SIZE*DISPLAY_SIZE) { map[i] = 0; /* blue */ } barrier(); /* Place each fish in a grid cell, then synchronize by waiting for all atomic operations to complete. */ for (i = 0, fish = fish_list; i < num_fish; i++, fish++) { x_disp = (int)(((fish->x_pos + DISPLAY_RANGE)*DISPLAY_SIZE)/ (2.0*DISPLAY_RANGE)); if (x_disp < 0 || x_disp >= DISPLAY_SIZE) continue; y_disp = (int)(((fish->y_pos + DISPLAY_RANGE)*DISPLAY_SIZE)/ (2.0*DISPLAY_RANGE)); if (y_disp < 0 || y_disp >= DISPLAY_SIZE) continue; /* SPLIT-C NOTE: The atomic routine performs an atomic operation on a piece of global data. The operation is performed on the processor to which the data is local. all_atomic_sync waits for all atomic operations to complete. */ atomic(add_fish, &map[x_disp + y_disp*DISPLAY_SIZE]); } all_atomic_sync(); return map; } /* Send changes in the ocean map between the argument and the last call's argument to the host. */ void all_display_fish(int *spread new_map) { static int init = 0; static int *spread old_map; int i; int x_disp, y_disp; /* If this is the first time, open a window, allocate the old map, and empty it. */ if (!init) { init = 1; on_one {X_init(DISPLAY_SIZE);} old_map = (int *spread) all_spread_malloc(DISPLAY_SIZE*DISPLAY_SIZE, sizeof(int)); /* SPLIT-C NOTE: The for_my_1d macro selects only those indices i for which the ith element of a spread array is local to the processor--this is equivalent to the 'owner computes' rule. */ for_my_1d(i, DISPLAY_SIZE*DISPLAY_SIZE) { old_map[i] = 0; } } /* Send changes in map to host for display and copy new map to old. */ for_my_1d(i, DISPLAY_SIZE*DISPLAY_SIZE) { if (new_map[i] != old_map[i]) { if (new_map[i] > 3) new_map[i] = 3; X_point(i%DISPLAY_SIZE, i/DISPLAY_SIZE, new_map[i]); } old_map[i] = new_map[i]; } } all_do_display(int num_fish, fish_t fish_list[]) { int *spread map; /* Generate display map. */ map = all_calculate_display(num_fish, fish_list); /* Send display data to host node. */ all_display_fish(map); /* Synchronize to insure that host has all of the current display data. */ barrier(); /* Display the ocean. */ on_one {X_show();} } /* Place fish in their initial positions. */ void all_init_fish(int num_fish, fish_t fish_list[]) { int i, n; double total_fish = PROCS*num_fish; fish_t *fish; for (i = 0, n = MYPROC*num_fish, fish = fish_list; i < num_fish; i++, n++, fish++) { fish->x_pos = n*2.0/total_fish - 1.0; fish->y_pos = 0.0; fish->x_vel = 0.0; fish->y_vel = fish->x_pos; fish->mass = 1.0 + n/total_fish; } } /* Add the force on fish1 due to fish2 to the force vector (fx, fy). */ void compute_force(fish_t *fish1, fish_t *fish2, double *fx, double *fy) { double x_sep, y_sep, dist_sq, grav_base; x_sep = fish2->x_pos - fish1->x_pos; y_sep = fish2->y_pos - fish1->y_pos; dist_sq = MAX((x_sep*x_sep) + (y_sep*y_sep), 0.01); /* Use the 2-dimensional gravity rule: F = d * (GMm/d^2) */ grav_base = GRAV_CONSTANT*(fish1->mass)*(fish2->mass)/dist_sq; (*fx) += grav_base*x_sep; (*fy) += grav_base*y_sep; } /* Move fish one time step. Update positions, velocity, and acceleration. Return local computations. */ void all_move_fish(int num_fish, fish_t *spread fish_list, double step, double *max_acc_ptr, double *max_speed_ptr, double *sum_speed_sq_ptr) { static double *spread x_force, *spread y_force; static int init = 0; int i, j; fish_t *fish; double *l_x_force, *l_y_force; double x_acc, y_acc; double cur_acc, max_acc = 0.0; double cur_speed, max_speed = 0.0; double speed_sq, sum_speed_sq = 0.0; /* Initialize spread arrays for force. */ if (!init) { init = 1; x_force = (double *spread)all_spread_malloc(NFISH, sizeof(double)); y_force = (double *spread)all_spread_malloc(NFISH, sizeof(double)); } /* Clear forces on local fish. */ for_my_1d(i, NFISH) { x_force[i] = 0.0; y_force[i] = 0.0; } /* First calculate force for local fish. */ for (j = 0; j < NFISH; j++) { fish_t rfish; rfish = fish_list[j]; for_my_1d(i, NFISH) { compute_force((fish_t *)&fish_list[i], &rfish, (double *)&x_force[i], (double *)&y_force[i]); } } /* Then move all local fish and return statistics. */ for (i = 0, fish = tolocal(fish_list), l_x_force = tolocal(x_force), l_y_force = tolocal(y_force); i < num_fish; i++, fish++, l_x_force++, l_y_force++) { /* Update fish positions, calculate acceleration, and update velocity. */ fish->x_pos += (fish->x_vel)*step; fish->y_pos += (fish->y_vel)*step; x_acc = (*l_x_force)/fish->mass; y_acc = (*l_y_force)/fish->mass; fish->x_vel += x_acc*step; fish->y_vel += y_acc*step; /* Accumulate local max speed, accel and contribution to mean square velocity. */ cur_acc = sqrt(x_acc*x_acc + y_acc*y_acc); max_acc = MAX(max_acc, cur_acc); speed_sq = (fish->x_vel)*(fish->x_vel) + (fish->y_vel)*(fish->y_vel); sum_speed_sq += speed_sq; cur_speed = sqrt(speed_sq); max_speed = MAX(max_speed, cur_speed); } /* Return local computation results. */ *max_acc_ptr = max_acc; *max_speed_ptr = max_speed; *sum_speed_sq_ptr = sum_speed_sq; } /* Simulate the movement of NFISH fish under a current. */ splitc_main() { int count = 1; double t = 0.0, dt = 0.01; double max_acc, max_speed, sum_speed_sq, mnsqvel; /* Allocate a global spread array for the fish data set. */ fish_t *spread fishes = all_spread_malloc(NFISH, sizeof(fish_t)); /* Obtain a pointer to the local portion of the spread array */ int num_fish = my_elements(NFISH); fish_t *fish_list = (fish_t *)&fishes[MYPROC]; all_init_fish(num_fish, fish_list); while (t < T_FINAL) { /* Display every STEPS_PER_DISPLAY time steps. */ if (--count == 0) { all_do_display(num_fish, fish_list); count = STEPS_PER_DISPLAY; } /* Update time. */ t += dt; /* Move fish with the current and compute rms velocity--notice that we send the global reference to the fish. */ all_move_fish(num_fish, fishes, dt, &max_acc, &max_speed, &sum_speed_sq); max_acc = all_reduce_to_all_dmax(max_acc); max_speed = all_reduce_to_all_dmax(max_speed); sum_speed_sq = all_reduce_to_all_dadd(sum_speed_sq); mnsqvel = sqrt(sum_speed_sq/NFISH); /* Adjust dt based on maximum speed and acceleration--this simple rule tries to insure that no velocity will change by more than 10% */ dt = 0.1*max_speed/max_acc; /* Print out time and rms velocity for this step. */ on_one {printf("%15.6lf %15.6lf;\n", t, mnsqvel);} } }