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