/* tab:4 * * wator.sc - particle-like fish under current * * Author: Steve Lumetta * Version: 1 * Creation Date: Wed May 26 15:34:46 1993 * Filename: wator.sc * History: * SL 1 Wed May 26 15:35:21 1993 * First written */ #include #include #include #include #include #define NFISH 10000 /* 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 */ /* Utility functions used in defining the current. */ #define DISTANCE(X,Y) (sqrt((X)*(X) + (Y)*(Y))) /* distance from origin */ #define MAX(X,Y) ((X) > (Y) ? (X) : (Y)) /* maximum of two numbers */ /* Two functions of x and y describe the external force due to the current. The functions here describe a whirlpool-like current. */ #define X_CURRENT(X,Y) (- 3.0*(Y)/MAX(DISTANCE(X,Y), 0.01) - X) #define Y_CURRENT(X,Y) (+ 3.0*(X)/MAX(DISTANCE(X,Y), 0.01) - Y) /* 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; } } /* Move fish one time step. Update positions, velocity, and acceleration. Return local computations. */ void all_move_fish(int num_fish, fish_t fish_list[], double step, double *max_acc_ptr, double *max_speed_ptr, double *sum_speed_sq_ptr) { int i; fish_t *fish; 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; for (i = 0, fish = fish_list; i < num_fish; i++, fish++) { /* Update fish positions, calculate acceleration, and update velocity. */ fish->x_pos += (fish->x_vel)*step; fish->y_pos += (fish->y_vel)*step; x_acc = X_CURRENT(fish->x_pos, fish->y_pos)/fish->mass; y_acc = Y_CURRENT(fish->x_pos, fish->y_pos)/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. */ all_move_fish(num_fish, fish_list, 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);} } }