/* 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 <split-c/split-c.h>
#include <split-c/control.h>
#include <split-c/com.h>
#include <split-c/atomic.h>
#include <math.h>
#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);}
}
}