/* tab:4
** fish2.c - particle-like fish under gravity
**
** Author: Xiaoye Li
** Version: 1
** Creation Date: Sat Jun 10 15:57:04 PDT 1995
** Filename: fish2.c
** History: 1. Origianl version did not have X output to display
** the moving fish. It appears to be bug free.
** 2. I added some routines to display the moving fish.
** The necessary files for linking were provided by the
** TA (Boris Vaysman - borisv@cs.berkeley.edu) and I
** did the actual modifications on fish1.c. As written,
** the fish are not erased when the new positions are
** plotted so fish leave a 'trail' as they move.
** Raymond Fellers (rsf@uclink3.berkeley.edu) 02/21/96
** 3. Fixed a race condition by inserting a barrier after
** all_init_fish (Boris Vaysman, borisv@cs)
**/
#include <thread.h>
#include <math.h>
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include "barrier.h"
#define NTHREADS 4 /* number of threads to do simulation */
#define NFISH 500 /* number of fish */
#define DISPLAY_SIZE 500 /* pixel size of display window */
#define SCALE 0.03 /* sets the magnification at the origin */
/* smaller #'s zoom in */
#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 0.01 /* proportionality constant of
gravitational interaction */
#define myID (thr_self()) /* my thread ID */
#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;
/*
* Information passed onto to a thread to tell it what to do.
*/
typedef struct
{
thread_t tid; /* thread ID */
int chunk; /* chunk of the global array assigned to this thread */
} thrinfo_t;
typedef struct {
volatile int zeroed; /* zeroed==1 means accum is already set to zero */
double accum;
} g_reduce_t;
/*
Globally shared variables
*/
fish_t *fishes;
double *x_force, *y_force;
thrinfo_t *thread_ptr = NULL;
mutex_t mul_lock;
barrier_t ba;
g_reduce_t g_dmax;
g_reduce_t g_dsum;
Display *theDisplay; /* These three variables are required to open the */
GC theGC; /* fish plotting window. They are externally */
Window theMain; /* declared in ui.h but are also required here. */
/*
Place fish in their initial positions.
*/
void all_init_fish(int mychunk, int num_fish, fish_t fish_list[])
{
int i, n;
double total_fish = NFISH;
fish_t *fish;
n = mychunk * num_fish;
fish = &fish_list[n];
for (i = 0; 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;
}
}
double all_reduce_to_all_dmax(double dmax)
{
barrier_wait(&ba);
if ( myID == thread_ptr[0].tid ) {
g_dmax.accum = 0.;
g_dmax.zeroed = 1;
}
while ( !g_dmax.zeroed ) ;
mutex_lock(&mul_lock);
g_dmax.accum = MAX(g_dmax.accum, dmax);
mutex_unlock(&mul_lock);
barrier_wait(&ba);
if ( myID == thread_ptr[0].tid )
g_dmax.zeroed = 0;
return (g_dmax.accum);
}
double all_reduce_to_all_dadd(double data)
{
barrier_wait(&ba);
if ( myID == thread_ptr[0].tid ) {
g_dsum.accum = 0.;
g_dsum.zeroed = 1;
}
while ( !g_dsum.zeroed ) ;
mutex_lock(&mul_lock);
g_dsum.accum += data;
mutex_unlock(&mul_lock);
barrier_wait(&ba);
if ( myID == thread_ptr[0].tid )
g_dsum.zeroed = 0;
return (g_dsum.accum);
}
/*
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 mychunk, int num_fish, fish_t *fish_list,
double step, double *max_acc_ptr,
double *max_speed_ptr, double *sum_speed_sq_ptr)
{
int i, j, k;
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;
/* Clear forces on local fish. */
k = mychunk * num_fish;
l_x_force = &x_force[k];
l_y_force = &y_force[k];
for (i = 0; i < num_fish; ++i) {
l_x_force[i] = 0.0;
l_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 (i = k; i < k + num_fish; ++i) {
compute_force((fish_t *)&fish_list[i], &rfish,
(double *)&x_force[i], (double *)&y_force[i]);
}
}
barrier_wait(&ba);
/* Then move all local fish and return statistics. */
for (i = 0, fish = &fish_list[k]; 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;
}
/*
Each thread begins execution by calling this function
*/
void *move_fish(void *arg_ptr)
{
int mychunk, num_fish;
int i, x, y;
fish_t *fish;
double t = 0.0, dt = 0.01;
double max_acc, max_speed, sum_speed_sq, mnsqvel;
mychunk = *(int*)arg_ptr;
num_fish = NFISH / NTHREADS;
all_init_fish(mychunk, num_fish, fishes);
barrier_wait(&ba); /* global synchronization is necessary to make
sure that positions of all fish have been
initialized */
while (t < T_FINAL) {
/* Update time. */
t += dt;
/* Move fish with the current and compute rms velocity. */
all_move_fish(mychunk, 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. */
if ( myID == thread_ptr[0].tid ) {
printf("%15.6lf %15.6lf;\n", t, mnsqvel);
}
/* Plot the movement of the fish */
/* clear_display (); uncomment this line if you don't want the fish*/
/* to leave a trail of dots. This is buggy... */
if (myID == thread_ptr[0].tid) {
fish = &fishes[0];
for (i=0; i<NFISH; i++, fish++) {
x = (int) (fish->x_pos/SCALE + DISPLAY_SIZE/2);
y = (int) (fish->y_pos/SCALE + DISPLAY_SIZE/2);
draw_point (x,y);
}
}
barrier_wait(&ba); /* don't let any procs. move the fish before */
/* thread 0 has finished plotting */
}
return 0;
}
/*
Simulate the movement of NFISH fish under a current.
*/
main()
{
int sync_type, i;
/* Allocate global shared arrays for the fish data set. */
fishes = (fish_t *) malloc(NFISH * sizeof(fish_t));
x_force = (double *) malloc(NFISH * sizeof(double));
y_force = (double *) malloc(NFISH * sizeof(double));
/* Initialize thread data structures */
thread_ptr = (thrinfo_t *) malloc(NTHREADS * sizeof(thrinfo_t));
sync_type = USYNC_PROCESS;
mutex_init(&mul_lock, sync_type, NULL);
barrier_init(&ba, NTHREADS, sync_type, NULL);
thread_ptr[0].chunk = 0;
thread_ptr[0].tid = myID;
for (i = 1; i < NTHREADS; i++) {
thread_ptr[i].chunk = i;
if (thr_create(0, 0, move_fish, (void*)&thread_ptr[i].chunk,
0, &thread_ptr[i].tid)) {
perror("thr_create");
exit(1);
}
}
/* Open an X window to display the fish */
simple_init (100,100,DISPLAY_SIZE,DISPLAY_SIZE);
/* Main thread starts simulation ... */
i = 0;
move_fish(&i);
/* Termination */
for (i = 1; i < NTHREADS; ++i)
thr_join(thread_ptr[i].tid, NULL, NULL);
printf("Hit return to close the window.");
getchar();
/* Close the X window used to display the fish */
XCloseDisplay(theDisplay);
return 0;
}