/* 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);}
}
}