/*                                  tab:4
** fish1.c - particle-like fish under current
**
** Author:          Xiaoye Li
** Version:         1
** Creation Date:   Sat Jun 10 13:14:19 PDT 1995
** Filename:        fish1.c
** History:         1. Original 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              100      /* number of fish */
#define DISPLAY_SIZE       500      /* pixel size of display window */
#define SCALE              0.1      /* 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 myID         (thr_self())   /* my thread ID */



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

/*
 * 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;
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 fishes[])
{
    int i, n;
    double total_fish = NFISH;
    fish_t *fish;

    n = mychunk * num_fish;
    fish = &fishes[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);
}

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

    fish = &fishes[mychunk * num_fish];
    for (i = 0; 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;
}

/*
   Each thread begins execution by calling this function
*/
void *move_fish(void *arg_ptr)
{
    int    mychunk, num_fish;
    int    count = 1;
    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. Buggy...            */

        /* allow only one processor to plot the fish */
        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;
    int    j, k;
    
    /* Allocate a global shared array for the fish data set. */
    fishes = (fish_t *) malloc(NFISH * sizeof(fish_t));

    /* 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 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 for fish plotting */
    XCloseDisplay(theDisplay);

    return 0;
}