nbody.c

一个在linux下的基于openmp的n-body编程

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<omp.h>
#define NUM_THREADS 5

typedef struct {
    double x;
    double y;
} coords;

struct body_s {
    coords p;
    coords v;
    coords f;
    double m;
    /* Body list pointers */
    struct body_s *next;
    struct body_s *prev;
};
typedef struct body_s body;

typedef struct {
  double px, py;
  double vx, vy;
  double fx, fy;
  double m;
} Body;


int my_rank;
int size;

Body* bodies;
int* myStart;
int* myEnd;
    

/* Global variables */
static int N;
static double width, height,G,DT,C;

/* List of bodies */
static body *blist_head = NULL;
static body *blist_tail = NULL;

/* Array with pointers to bodies */
static body **barray;

/* Insert a node with value v into the body list.
 */
int blist_insert(body *new)
{
    /*
     * Insert node at the tail of the list.
     */
    new->next = NULL;
    new->prev = blist_tail;
    if (blist_tail != NULL)
	blist_tail->next = new;
    if (blist_head == NULL) 
	blist_head = new;
    blist_tail = new;
	
    return 0;
}

int readline(FILE *fp, void *vptr, size_t maxlen) {
    int n, rc;
    char c, *ptr;

    ptr = vptr;
    /* n = 1, because we need space to store '\0' */
    for (n = 1; n < maxlen; n++) {
	if ((rc = fread(&c, 1, 1, fp)) == 1) { /* one byte read */
	    *ptr++ = c;

	    if (c == '\n') /* end of line */
		break;
	}
	else if (rc == 0) {
	    if (n == 1)
		/* EOF, no data read */
		return 0;
	    else
		/* EOF, data read */
		break;
	}
    }
    
    *ptr = '\0';
    return n;
}

#define LINE_SIZE 250   /* Maximum size of a line in the data file */
#define BODY_ENTRIES 7  /* Number of entries per body (as described belov) */

void read_input(const char *filename) {
    FILE *fp = fopen(filename, "r");
    char line[LINE_SIZE], *p, *p_new;
    int current_line = 0;
    body *new;
    int i;
    body *bp;
    
    /******************************************/
    /*     Read base information              */
    /******************************************/
    
    /*read N*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%d", &N);
	printf("N:%d\n",N);
    /*read G*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%lf", &G);
	printf("G:%lf\n",G);
    /*read Size*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%lf %lf", &width, &height);
	printf("Size:%lf %lf\n",width,height);
    /*read DT*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%lf", &DT);
	printf("DT:%lf\n",DT);
    /*read c*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%lf", &C);
	printf("C:%lf\n",C);
    /* Number of bodies */

    
    while (readline(fp, (void *)line, LINE_SIZE) != 0) {

	/* Check line for correctnes (number of entries) */
	current_line++;
//	if (check_line(line) == -1)
//	    printf("Error at line %d: %s\n", current_line, line);

	new = malloc(sizeof(body));

	/* Parse line */
	p = line;
	/* Goto next non-blank character */
	while (isspace(*p) && *p != '\0') {
	    *p = '\0';
	    p++;
	    printf(".");
	}

	/* Set entries */
	new->m = strtod(p, &p_new);
	p = p_new;
	new->p.x = strtod(p, &p_new);
	p = p_new;
	new->p.y = strtod(p, &p_new);
	p = p_new;
	new->v.x = strtod(p, &p_new);
	p = p_new;
	new->v.y = strtod(p, &p_new);
		
    
	/* Insert in list */
	blist_insert(new);
    }	
    
    //N = current_line; /* Number of bodies */
    bodies = (Body *) malloc(sizeof(Body) * N);
    for (i = 0, bp = blist_head; i < N; i++, bp = bp->next) {
      bodies[i].px = bp->p.x;
      bodies[i].py = bp->p.y;
      bodies[i].vx = bp->v.x;
      bodies[i].vy = bp->v.y;
      bodies[i].fx = 0.0;
      bodies[i].fy = 0.0;
      bodies[i].m = bp->m;
    }

    /* Initialize barray */
    barray = malloc(sizeof(body *) * N);
    for (i = 0, bp = blist_head; i < N; i++, bp = bp->next)
	barray[i] = bp;
}

/* Calculate forces between my part of the bodies */
void calculate_forces(void) {
  int i, j;
  coords direction;
  double distance,distance2, magnitude;
  omp_set_num_threads(NUM_THREADS);
#pragma omp parallel for default(shared) private(i,j,direction,distance,distance2,magnitude)
  for(i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      if (i == j) continue;
	  distance2 = pow((bodies[i].px - bodies[j].px), 2.0) + 
		       pow((bodies[i].py - bodies[j].py), 2.0);
      distance = sqrt( distance2 );
      magnitude = (G * bodies[i].m * bodies[j].m) / (distance2+pow(C, 2.0));
      direction.x = bodies[j].px - bodies[i].px;
      direction.y = bodies[j].py - bodies[i].py;

      bodies[i].fx = bodies[i].fx + magnitude * direction.x / distance;
      bodies[i].fy = bodies[i].fy + magnitude * direction.y / distance;
    }
  }
}

void move_bodies(void) {
  int i;
  coords deltav;
  coords deltap;
  double tempPointx,tempPointy;
  omp_set_num_threads(NUM_THREADS);
  #pragma omp parallel for default(shared) private(i,deltav,deltap,tempPointx,tempPointy)
  for (i = 0; i < N; i++) {
    deltav.x = bodies[i].fx / bodies[i].m * DT;
    deltav.y = bodies[i].fy / bodies[i].m * DT;
    deltap.x = (bodies[i].vx + deltav.x / 2) * DT;
    deltap.y = (bodies[i].vy + deltav.y / 2) * DT;

    tempPointx = bodies[i].px + deltap.x;
    tempPointy = bodies[i].py + deltap.y;
	//运行轨迹超出了空间范围
	if(tempPointx < 0 || tempPointx > width ||
	tempPointy < 0 || tempPointy > height){
	bodies[i].vx = -bodies[i].vx;
    bodies[i].vy = -bodies[i].vy;
    bodies[i].px = bodies[i].px;
    bodies[i].py = bodies[i].py;
	}
	else
	{
    bodies[i].vx = bodies[i].vx + deltav.x;
    bodies[i].vy = bodies[i].vy + deltav.y;
    bodies[i].px = tempPointx;
    bodies[i].py = tempPointy;
	}

    /* Reset force vectors */
    bodies[i].fx = 0.0;
    bodies[i].fy = 0.0;

  }
}

/* Run simulation for 'steps' steps. Results are writen to res_filename. */
void run(int steps, char *res_filename) {
  int step, j,i;
  FILE *fp;
    printf("Running %i bodies, %i steps\n", N, steps);
    fp = fopen(res_filename, "w");

    /* Calculate */
	for(i=0;i<steps;i++)
	{
	printf("steps:%d\n",i);
    calculate_forces();
    move_bodies();
	}


//输出结果;
      fprintf(fp, "%d\n",N);
      fprintf(fp, "%lf\n",G);
      fprintf(fp, "%lf %lf\n",width,height);
      fprintf(fp, "%lf\n",DT);
      fprintf(fp, "%lf\n",C);
      for (j = 0; j < N; j++)
      fprintf(fp, "%lf %lf %lf %lf %lf\n",bodies[j].m, bodies[j].px, bodies[j].py, bodies[j].vx, bodies[j].vy);
    fclose(fp);
  
}

int main(int argc, char *argv[]) {
  int i, z;
  read_input("sample_input.in");
  run(1000, "sample_out.out");
  return 0;
}