///////////////////////////////////////////////////////////////////////////////
// $Id: fish.c,v 1.5 1998/02/20 19:25:29 dmartin Exp $
///////////////////////////////////////////////////////////////////////////////

#include <math.h>
#include <stdlib.h>
#include "dmartin/tablet.h"
#include "dmartin/timer.h"
#include "mpi.h"

///////////////////////////////////////////////////////////////////////
// MACROS
///////////////////////////////////////////////////////////////////////

#define MAX(X,Y) ((X) > (Y) ? (X) : (Y))  
#define MIN(X,Y) ((X) > (Y) ? (Y) : (X))  
#define SQR(X)   ((X)*(X))
#define ABS(X)   ((X>=0) ? X : -(X))
#define DISTANCE(X,Y)  (sqrt((X)*(X) + (Y)*(Y)))

///////////////////////////////////////////////////////////////////////
// TYPES
///////////////////////////////////////////////////////////////////////

typedef enum {true=1, false=0} bool;

typedef struct {
  double d_px, d_py;  // position
  double d_vx, d_vy;  // velocity
  double d_m;         // mass
} fish_t;

///////////////////////////////////////////////////////////////////////
// GLOBALS
///////////////////////////////////////////////////////////////////////

static int g_i_id;

///////////////////////////////////////////////////////////////////////
// ARGUMENT PROCESSING
///////////////////////////////////////////////////////////////////////

static void
usage (const char* pc_prog)
{
  if (g_i_id == 0)
    {
      fprintf (stderr, 
	       "usage: %s "
	       "-fish <n> "
	       "-procs <n> "
	       "-dtmin <x> "
	       "-dtmax <x> "
	       "-dxmin <x> "
	       "-grav <x> "
	       "-end <x> "
	       "-display <yes|no> "
	       "-size <n>"
	       "\n", pc_prog);
    }
  exit (1);
}

static void
badParam (const char* pc_prog)
{
  if (g_i_id == 0)
    {
      fprintf (stderr, "%s: bad parameter value\n", pc_prog);
    }
  exit (1);
}

static int 
intArg (char* pc_prog, int* argc, char*** argv)
{
  int ret;
  (*argc)--; (*argv)++;
  if (*argc == 0) usage (pc_prog);
  ret = atoi (**argv);
  (*argc)--; (*argv)++;
  return ret;
}

static double
floatArg (char* pc_prog, int* argc, char*** argv)
{
  double ret;
  (*argc)--; (*argv)++;
  if (*argc == 0) usage (pc_prog);
  ret = atof (**argv);
  (*argc)--; (*argv)++;
  return ret;
}

static bool
boolArg (char* pc_prog, int* argc, char*** argv)
{
  bool ret;
  (*argc)--; (*argv)++;
  if (*argc == 0) usage (pc_prog);
  if (strcmp (**argv, "yes") == 0)
    {
      ret = true;
    }
  else if (strcmp (**argv, "no") == 0)
    {
      ret = false;
    }
  else
    {
      usage (pc_prog);
    }
  (*argc)--; (*argv)++;
  return ret;
}

///////////////////////////////////////////////////////////////////////////////
// MAIN
///////////////////////////////////////////////////////////////////////////////

int
main (int argc, char **argv)
{
  MPI_Datatype MPI_fish;
  MPI_Status status;
  int i_procs, i_fish, i_windowSize;
  double d_t, d_dt, d_dtMin, d_dtMax, d_endTime, d_gravConst, d_dxMin;
  double d_sqrSpeed, d_maxSqrSpeed, d_myMaxSqrSpeed;
  bool b_display;
  char* pc_prog = *argv;
  int i_tag = 50;
  int i, j;
  double* pd_fx;
  double* pd_fy;
  int* pi_fishStart;
  int* pi_fishEnd;
  fish_t* p_fish;
  void* p_window;
  void* p_timer;

  // pass arguments to MPI
  MPI_Init (&argc, &argv);

  // get my processor id
  MPI_Comm_rank (MPI_COMM_WORLD, &g_i_id);          

  // get number of processors
  MPI_Comm_size (MPI_COMM_WORLD, &i_procs);       

  // argument defaults
  i_fish = 50;
  d_endTime = 50;
  d_dtMin = 0.0001;
  d_dtMax = 0.01;
  d_dxMin = 0.005;
  d_gravConst = 1.0;
  b_display = true;
  i_windowSize = 200;
  
  // process arguments
  argc--; argv++;
  while (argc > 0)
    {
      if (strcmp (*argv, "-fish") == 0) 
	{
	  i_fish = intArg (pc_prog, &argc, &argv);
	}
      else if (strcmp (*argv, "-display") == 0) 
	{
	  b_display = boolArg (pc_prog, &argc, &argv);
	}
      else if (strcmp (*argv, "-dtmin") == 0) 
	{
	  d_dtMin = floatArg (pc_prog, &argc, &argv);
	}
      else if (strcmp (*argv, "-dtmax") == 0) 
	{
	  d_dtMax = floatArg (pc_prog, &argc, &argv);
	}
      else if (strcmp (*argv, "-dxmin") == 0) 
	{
	  d_dxMin = floatArg (pc_prog, &argc, &argv);
	}
      else if (strcmp (*argv, "-grav") == 0) 
	{
	  d_gravConst = floatArg (pc_prog, &argc, &argv);
	}
      else if (strcmp (*argv, "-end") == 0) 
	{
	  d_endTime = floatArg (pc_prog, &argc, &argv);
	}
      else if (strcmp (*argv, "-size") == 0) 
	{
	  i_windowSize = intArg (pc_prog, &argc, &argv);
	}
      else
	{ 
	  usage (pc_prog);
	}
    }

  // check argument values
  if (i_fish < 1 || 
      i_procs > i_fish ||
      d_dtMin > d_dtMax ||
      d_dxMin < 0 ||
      d_endTime < 0 ||
      i_windowSize < 1) 
    {
      badParam (pc_prog);
    }

  // report simulation parameters
  if (g_i_id == 0)
    {
      fprintf (stderr, "procs = %d\n", i_procs);
      fprintf (stderr, "fish = %d\n", i_fish);
      fprintf (stderr, "min dt = %f\n", d_dtMin);
      fprintf (stderr, "max dt = %f\n", d_dtMax);
      fprintf (stderr, "min dx = %f\n", d_dxMin);
      fprintf (stderr, "gravitational constant = %f\n", d_gravConst);
      fprintf (stderr, "simulation end = %f\n", d_endTime);
      fprintf (stderr, "display = %s\n", b_display ? "yes" : "no");
      fprintf (stderr, "window size = %d\n", i_windowSize);
    }

  // let MPI know about fish_t
  MPI_Type_contiguous (5, MPI_DOUBLE, &MPI_fish);
  MPI_Type_commit (&MPI_fish);

  // allocate stuff
  p_fish = (fish_t *) malloc (i_fish * sizeof (fish_t));
  pi_fishStart = (int *) malloc (i_procs * sizeof (int));
  pi_fishEnd = (int *) malloc (i_procs * sizeof (int));
  pd_fx = (double*) malloc (i_fish * (sizeof (double)));
  pd_fy = (double*) malloc (i_fish * (sizeof (double)));
  if (b_display && g_i_id == 0)
    {
      p_window = Ctablet_new ();
      Ctablet_init (p_window, "Fish!", i_windowSize, i_windowSize);
    }
  p_timer = Ctimer_new ();
  
  // calculate which fish belong to which processor
  pi_fishStart[0] = 0;
  pi_fishEnd[0] = i_fish / i_procs;
  for (i = 1; i < i_procs; i++)
    {
      pi_fishStart[i] = pi_fishEnd[i-1];
      pi_fishEnd[i] = pi_fishStart[i] + i_fish / i_procs;
    }
  pi_fishEnd[i_procs-1] = i_fish; 
  if (g_i_id == 0)
    {
      for (i = 0; i < i_procs; i++)
	{
	  //fprintf (stderr, "processor %d owns fish [%d..%d)\n", i, pi_fishStart[i], pi_fishEnd[i]);
	}
    }

  // initialize my fish 
  for (i = pi_fishStart[g_i_id]; i < pi_fishEnd[g_i_id]; i++) 
    {
      double z = 2.0 * M_PI / i_fish;

      p_fish[i].d_px = sin (z * i);
      p_fish[i].d_py = cos (z * i);
      p_fish[i].d_vx = 0;
      p_fish[i].d_vy = 0;
      p_fish[i].d_m  = 1.0 / i_fish;
    }

  // start timer
  if (g_i_id == 0)
    {
      Ctimer_start (p_timer);
    }

  // loop over the time steps 
  for (d_dt = d_dtMin, d_t = 0; d_t < d_endTime; d_t += d_dt)
    {
      // all-to-all exchange of fish information
      // ...first gather everything to processor 0 
      if (g_i_id == 0) 
	{
	  for (i = 1; i < i_procs; i++) 
	    {
	      MPI_Recv (&p_fish[pi_fishStart[i]], 
			pi_fishEnd[i] - pi_fishStart[i],
			MPI_fish, i, i_tag, MPI_COMM_WORLD, &status); 
	    }
	}
      else 
	{
	  MPI_Send (&p_fish[pi_fishStart[g_i_id]], 
		    pi_fishEnd[g_i_id] - pi_fishStart[g_i_id],
		    MPI_fish, 0, i_tag, MPI_COMM_WORLD);
	}

    // ...now broadcast the fish array
    MPI_Bcast (p_fish, i_fish, MPI_fish, 0, MPI_COMM_WORLD);

    // compute forces on my fish 
    for (i = pi_fishStart[g_i_id]; i < pi_fishEnd[g_i_id]; i++) 
      {
	// compute global force on fish i
	pd_fx [i] = 0;
	pd_fy [i] = 0;
	for (j = 0; j < i_fish; j++) 
	  {
	    double d_dx = p_fish[j].d_px - p_fish[i].d_px;
	    double d_dy = p_fish[j].d_py - p_fish[i].d_py;
	    double d_sqrDist = MAX (SQR(d_dx) + SQR(d_dy), SQR(d_dxMin));
	    double d_gravBase = d_gravConst * p_fish[i].d_m * p_fish[j].d_m / d_sqrDist;
	    pd_fx[i] += d_dx * d_gravBase;
	    pd_fy[i] += d_dy * d_gravBase;
	  }
      }

    // compute new positions and velocities of my fish
    // compute maximum velocity of my fish
    d_myMaxSqrSpeed = 0;
    for (i = pi_fishStart[g_i_id]; i < pi_fishEnd[g_i_id]; i++) 
      {
	p_fish[i].d_px += p_fish[i].d_vx * d_dt;
	p_fish[i].d_py += p_fish[i].d_vy * d_dt;
	p_fish[i].d_vx += pd_fx[i] / p_fish[i].d_m * d_dt;
	p_fish[i].d_vy += pd_fy[i] / p_fish[i].d_m * d_dt;
	d_sqrSpeed = SQR (p_fish[i].d_vx) + SQR (p_fish[i].d_vy);
	d_myMaxSqrSpeed = MAX (d_myMaxSqrSpeed, d_sqrSpeed);
      }  

    // compute global maximum fish velocity 
    MPI_Allreduce (&d_myMaxSqrSpeed, &d_maxSqrSpeed, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
    
    // recompute dt so that the fastest particle does not move more
    // than the minumum distance in the next time step
    d_dt = d_dxMin / sqrt (d_maxSqrSpeed);
    d_dt = MIN (d_dtMax, d_dt);
    d_dt = MAX (d_dtMin, d_dt);
    
    // display fish
    if (b_display && g_i_id == 0)
      {
	for (i = 0; i < i_fish; i++)
	  {
	    Ctablet_drawPoint (p_window, 
			       (p_fish[i].d_px + 1) / 2 * 0.5 + 0.25,
			       (p_fish[i].d_py + 1) / 2 * 0.5 + 0.25);
	  }
	Ctablet_renderScene (p_window);
      }
    }

  // stop timer and display time
  if (g_i_id == 0)
    {
      Ctimer_stop (p_timer);
      fprintf (stderr, "time = %f sec\n", Ctimer_value (p_timer));
    }

  MPI_Finalize();
  return 0;
}

///////////////////////////////////////////////////////////////////////////////
// $Log: fish.c,v $
// Revision 1.5  1998/02/20 19:25:29  dmartin
// seems to work
//
// Revision 1.4  1998/02/09 20:10:31  dmartin
// *** empty log message ***
//
// Revision 1.3  1998/02/09 19:48:00  dmartin
// *** empty log message ***
//
// Revision 1.2  1998/02/09 19:45:05  dmartin
// *** empty log message ***
//
// Revision 1.1  1998/02/09 19:41:53  dmartin
// Initial revision
//
///////////////////////////////////////////////////////////////////////////////