/*                   VERSION 4.0
		     05/13/2002

   This version has 3 different kind of fields
   implemented and all the bugs in odestep.for,
   emflds?.for has been fixed.

   This version inlcudes a new common block /times
   in loadcommon.for.

   The setup file was modified to include dt as an
   input in the setup file.

   A new file progress_file is updated on node 0 to
   indicate what the program is working on, thus 
   enabling the users to check what the program is working
   on at any given instant of time

   New subroutine error_printf__(char *display_string)
   has been added and this is used by odestep.for and 
   rkcs.for to denote any Step Size Underflow.

		     Version 3.0
   This is the version that implements different method
   for reading in the varies parameter making use of
   a new subroutine called loadcommon.for

                     Version 2.1
   The newest addition is the type of derivative and 
   the type of field could be individually selected
   this is done by setting up the following two variables
   in the setup file

   derivative_type = 1 or 2 
        1 -> field lines
	2 -> Charge particles

   and field_type = 1 or 2.
	1 -> 1d code
	2 -> 2d code
	3 -> 
*/

/* This program still carries the previous verison's option of
   writing the output file in node 0 it also include the 
   original option of writing the output file in
   individual nodes option changed by the following variable
   in the setup file:

   output_node_0 = 0 or 1
        0 -> write it in individual nodes
	1 -> write it in node 0 in a single file
                                                 */

/* program Streamline	
   Pushes  charged test particles in prescribed EM fields
   This version adaptively adjusts the timestep to maintain
   accuracy, in accordance with the particle with the worst
   error for each timestep. The turbulent fields are 
   artificially generated according to a specified spectral
   distribution. No nonlinear wave-wave interactions are
   accounted for. This version uses 1d grids for the EM 
   fields.                    */

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>

#define DEBUG  0
#define DEBUG1 0
#define DEBUG2 0
#define ILIMIT 50

/*  NX, NY, NZ are the number of grid points in the x, y 
    and z directions, for the electromagnetic field arrays */

int NZ, NX, NY;

/* Various environment variable */
int derivative_type1, field_type1, nfs, output_node_0;

/* Following are the new modifications suggested */

double *job_input_array, *job_result;

char output_filename[80];

double eps, racc, Det_t, dbob1, b01;
double lenx1, leny1, lenz1, dt;
int nits, num_fieldlines;

int num_integers, num_doubles;

int temp_int_transfer[8];
double temp_double_transfer[8];

int integer_input[ILIMIT];
double double_input[ILIMIT];

/* int tag; */

int length = 10;       /* the length of the job array */

main( int argc, char **argv)
{
  /* Variables for MPI use */
  int my_rank, p, s, bg = 0, np, current_job = 0, completed;

  int go_on, root_loop; /* variable for child status check */
  extern int length;    /* length of the job array */

  FILE *setup_file, *output_file, *assign_record;
 
  double start_time, end_time;
  double *job;   /*  Input and output data arrays */

  MPI_Status status;
  int fnsh=0, fnsh_proc=0;
  
  char *out_file;

  int i;

  /* Starting of MPI Initialization */
  MPI_Init( &argc, &argv);
  
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  
  MPI_Comm_size(MPI_COMM_WORLD, &np);

  start_time=MPI_Wtime();

  out_file = (char *) malloc (80 * sizeof(char) );

  /* End of MPI Initialization */

  if ( ( setup_file = fopen("setup", "r") ) == NULL )
    {
      printf("File not found: setup\nProgram Aborted\n");
      MPI_Abort(MPI_COMM_WORLD,0);
      exit(0);
    }

    fscanf(setup_file, "NFS=%d\n", &nfs);
  fclose(setup_file);     

  if(nfs) /* If NFS is present every node reads the file .... more efficient*/
    input_file_read( my_rank);

  else
    {
      /* If NFS not present node 0 reads the file and broadcast to everynode */

      if ( my_rank == 0 )
	input_file_read( my_rank);      

      /* Assignment of values to this temp variable 
	 is done within input_file_read() */
      
      MPI_Bcast(temp_int_transfer, 8, MPI_INT, 0, MPI_COMM_WORLD);
      
      MPI_Bcast(temp_double_transfer, 8, MPI_DOUBLE, 0, MPI_COMM_WORLD);

      MPI_Bcast(integer_input, ILIMIT, MPI_INT, 0, MPI_COMM_WORLD);

      MPI_Bcast(double_input, ILIMIT, MPI_DOUBLE, 0, MPI_COMM_WORLD);

      MPI_Bcast(output_filename, 80, MPI_CHAR, 0, MPI_COMM_WORLD);

      eps = racc;

      if ( my_rank != 0)
	{
	  Det_t  = temp_double_transfer[0];
	  racc   = temp_double_transfer[1];
	  dbob1  = temp_double_transfer[2];
	  b01    = temp_double_transfer[3];
	  lenx1  = temp_double_transfer[4];
 	  leny1  = temp_double_transfer[5];
   	  lenz1  = temp_double_transfer[6];
	  dt     = temp_double_transfer[7];
   
	  output_node_0  = temp_int_transfer[0];
	  derivative_type1= temp_int_transfer[1];
          field_type1    = temp_int_transfer[2];
	  num_fieldlines = temp_int_transfer[3];
	  NZ             = temp_int_transfer[4];
	  NX             = temp_int_transfer[5];
	  NY             = temp_int_transfer[6];
	  nits           = temp_int_transfer[7];
	}
    }

  /* Commands to be executed by the master */
  if ( my_rank == 0 )
    {
      /* Memory allocation for the job array */
      if ( (job_input_array = (double *) malloc( length * num_fieldlines * sizeof(double) ))== NULL )
	{
	  printf("Failed memory allocation: job_input_array\nProgram Aborted\n");
	  MPI_Abort(MPI_COMM_WORLD,0);
	  exit(0);
	}
      
      /* jobs are read from the file data and 
	      job array is intialized in here */
      x0_file_read();
    }

  loadcommon_( &NX, &NY, &NZ, &lenx1, &leny1, &lenz1, \
	       &Det_t, &nits, &dt, \
	       &derivative_type1, &field_type1, \
	       &racc, &dbob1, &b01, integer_input, double_input, \
	       &num_integers, &num_doubles);

  sprintf(out_file, "%s%d", output_filename, my_rank);

  /* setting up of FILES */
  if (output_node_0)
    {
      /*Write in individual node with filename = output_filename+my_rank */
      if (my_rank == 0)
	if ( ( output_file = fopen (out_file, "w") ) == NULL)
	  {
	    printf("Unable to create file: out_file\nProgram Aborted\n");
	    MPI_Abort(MPI_COMM_WORLD,0);
	    exit(0);
	  }
    }
  else
    {
      /* Write in node 0 with filename = output_filename */
      if (my_rank != 0)
	if ( ( output_file = fopen (out_file, "w") ) == NULL)
	  {
	    printf("Unable to create file: out_file\nProgram Aborted\n");
	    MPI_Abort(MPI_COMM_WORLD,0);
	    exit(0);
	  }
    }

  if ( my_rank != 0 ) /* code for all the number crunching nodes*/
    {  
      go_on = 1;
      bg = my_rank;
      
      while ( go_on ) /* Infinite loop till go_on is set to 0 */
	{
	  if (DEBUG) printf("from WHILE loop of node %d\n", my_rank);
	  
	  /* intialize the vaarious arrays */
	  job = (double *)malloc( length * sizeof(double));
	  job_result= (double *)malloc( (nits + 1) * length * sizeof(double));

	  if (job == NULL || job_result == NULL)
	    {
	      printf("Failed memory allocation: job and job_result\nProgram Aborted\n");
	      MPI_Abort(MPI_COMM_WORLD,0);
	      exit(0);
	    }
	  /* Inform node 0 ..... the availability */
	  MPI_Send( &bg, 1, MPI_INTEGER, 0, my_rank, MPI_COMM_WORLD);

	  if (DEBUG) printf("Past the send statement\n");

	  /* Recieve from the node 0 .... the next job or end condition */
	  MPI_Recv( &job[0], length, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status );
	  
	  if (DEBUG) printf("Past the send and recv statements in node %d\n", my_rank);

	  if ( job[0] > 0.0 ) /* Check for the end condition */
	    /* The end codition is job[0] == -1 */
	    {
	      /* call myodestep which computes the whole field line */
	      completed = myodestep( job, job_result);

	      /* completed gives the no. of iterations completed before 
		 any error occured this is used with local_file_write
		 in case of local_file_write ... in node 0 writing 
		 nits -> no of entries are written as default
		 completed is not used in local_file_write */

	      if (DEBUG) printf("Done with myodestep");

	      if (output_node_0)

		/* Send a message to node 0 ---- with the job_result array*/
		MPI_Send(job_result, (nits + 1)*length, MPI_DOUBLE, 0, my_rank, MPI_COMM_WORLD);

	      else
		/* Write the output to a local file */
		local_file_write( completed, job_result, output_file);

	      if (DEBUG) printf("Done with file writing\n");

	      free(job);
	      free(job_result);
	    }

	  else
	    {
	      go_on = 0;
	      /* printf("ELSE part\n"); */
	    }
	}
      if (DEBUG) printf("\nChild Node %d exiting\n", my_rank);
    }

  /* End of other nodes code start of the Master's code section */
  else
    {
      
      /* Deleting the old progress_file and creating a new file */
      if ( ( assign_record = fopen ("progress_file", "w") ) == NULL)
	{
	  printf("Unable to create file : progress_file\nProgram Aborted\n");
	  MPI_Abort(MPI_COMM_WORLD,0);
	  exit(0);
	}
      
      fclose (assign_record);

      for( root_loop = 0; root_loop < np - 1; root_loop++) 
	/* removed notfinished routine here with the whole function */
	/* Dispatch of the first set of job to each node */
	{
	  if ( (job = (double *)malloc( length * sizeof(double)) ) == NULL )
	    {
	      printf("Failed memory allocation: job\nProgram Aborted\n");
	      MPI_Abort(MPI_COMM_WORLD,0);
	      exit(0);
	    }

	  /* Allocate memory for the job array and get the next job */

	  getnewjob(job, job_input_array, num_fieldlines, current_job);
	  current_job++; /* counter for no. of jobs */

	  if (DEBUG1) printf("1 Before Recv\n");

	  MPI_Recv( &bg, 1, MPI_INTEGER, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
	  
	  if (DEBUG1) printf("1 past recv and before send in node 0\n");	  
	  
	  /* bg contains the rank of the node that sends out the request  */	  	  

	  MPI_Send( &job[0], length, MPI_DOUBLE, bg, 0, MPI_COMM_WORLD);

	  if (job[0] > 0.0)
	    {
	      assign_record = fopen ("progress_file", "a") ;
	      fprintf( assign_record, "Job Number %d has been assigned to node %d\n", current_job, bg);
	      fclose(assign_record);
	    }

	  if (DEBUG1) printf("1 Past the first send\n");
	  free(job);
	}

      if (DEBUG1) printf("Flag %d %d\n", np, num_fieldlines + np - 1);

      for( root_loop = np-1; root_loop < num_fieldlines + np - 1; root_loop++) 
	/* Rest of the job is issued from within this loop */
	{
	  job_result= (double *)malloc( (nits + 1) * length * sizeof(double));
	  
	  job = (double *)malloc( length * sizeof(double));

	  if (job == NULL || job_result == NULL)
	    {
	      printf("Failed memory allocation: job and job_result\nProgram Aborted\n");
	      MPI_Abort(MPI_COMM_WORLD,0);
	      exit(0);
	    }

	  getnewjob(job, job_input_array, num_fieldlines, current_job);
	  current_job++;

	  if (DEBUG2) printf("\njob[0] = %f\n", job[0]);

	  /* Dispatch of the first set of job to each node */
	  
	  if (output_node_0)
	    {
	      /* Recieve the result from the other nodes */

	      /* There is a limitation here nits*length should be within
		 the legal limits of the integer else the count becomes -ve
		 and it seg faults */

	      MPI_Recv( job_result, (nits+1) * length, MPI_DOUBLE, 
			MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
	      
	      local_file_write( (nits + 1), job_result, output_file);

	      /* Recieve the next job request from the same node from  which
		 job_result was recieved together with the same tag */
	      MPI_Recv( &bg, 1, MPI_INTEGER, status.MPI_SOURCE, status.MPI_TAG,
		    MPI_COMM_WORLD, &status);
	    }
 
	  else
	    /* Recieve the request for the next job */
	    MPI_Recv( &bg, 1, MPI_INTEGER, MPI_ANY_SOURCE, MPI_ANY_TAG,
		      MPI_COMM_WORLD, &status);
	  
	  /* bg contains the rank of the node that sends out the request  */	  	  
	  
	  MPI_Send( &job[0], length, MPI_DOUBLE, bg, 0, MPI_COMM_WORLD);

	  if (job[0] > 0.0)
	    {
	      assign_record = fopen ("progress_file", "a") ;
	      fprintf( assign_record, "Job completed on node %d\n", bg);
	      fprintf( assign_record, "Job Number %d has been assigned to node %d\n", current_job, bg);
	      fclose(assign_record);
	    }
	  
	  if (DEBUG) printf("2 past send in node 0\n");
	 
	  free(job);
	  free(job_result);
	}

      if (DEBUG) printf("\nNode 0 exiting\n");
    }
  
  end_time = MPI_Wtime();
  if (my_rank == 0)
    printf("# of fldlines = %d\ntotal time = %f\n# of processors = %d\n", 
	   num_fieldlines, end_time - start_time, np); 
  
  free(job_input_array);
  if (output_node_0)
    {
      if (my_rank == 0)
	fclose(output_file);
    }
  else
    if ( my_rank != 0 )
      fclose(output_file);
  MPI_Finalize();
}







