First working commit

Makefile

+all: hpprp run density/video.avi
+
+hpprp: hpprp.c config.c
+	g++ -O3 hpprp.c -o hpprp
+
+run: hpprp
+	./hpprp
+
+density/video.avi: run
+	./image_gen.sh density
+
+clean:
+	rm -f density/*.txt density/*.png density/video.avi hpprp

config.c

+// USER CONFIGURATION FILE
+
+// Main properties
+#define WIDTH 400
+#define HEIGHT 300
+#define ITERS 1000
+#define SAVE_PERIOD 10
+#define AVERAGING_RADIUS 5
+
+// Field initialization
+void init()
+{
+	int i, j;
+	for (i=0; i<HEIGHT; i++) {
+		for (j=0; j<WIDTH; j++) {
+			double r=drand48();
+			if (r>=0 && r<1.0/3) {
+				FIELD(i, j)=5;
+			} else if (r>=1.0/3 && r<2.0/3) {
+				FIELD(i, j)=10;
+			} else {
+				FIELD(i, j)=16;
+			}
+		}
+	}
+}
+
+// Collision rules
+char collide(char cell)
+{
+	double r=drand48();
+	switch(cell) {
+		case 5:
+		case 10:
+		case 16:
+			if (r>=0 && r<1.0/3) { return 5; }
+			else if (r>=1.0/3 && r<2.0/3) { return 10; }
+			else { return 16; }
+		default:
+			return cell;
+	}
+}
+
+// Sources rules
+void sources(int iter)
+{
+	if (iter==0) {
+		int i, j;
+		for (i=0; i<HEIGHT; i++) {
+			for (j=WIDTH*4/10; j<WIDTH*6/10; j++) {
+				FIELD(i, j)=15;
+			}
+		}
+	}
+}
+

hpprp.c

+#include <omp.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#define FIELD(i, j) field[(i)*(WIDTH) + j]
+#define NEW_FIELD(i, j) new_field[(i)*(WIDTH) + j]
+char *field, *new_field;
+
+#include "config.c"
+inline int get_mass(char cell)
+{
+	int mass=0;
+	if (cell&1) mass+=1;
+	if (cell&2) mass+=1;
+	if (cell&4) mass+=1;
+	if (cell&8) mass+=1;
+	if (cell&16) mass+=2;
+	if (cell&32) mass+=4;
+	if (cell&64) mass+=8;
+	if (cell&128) mass+=16;
+	return mass;
+}
+
+int sum_mass(int row, int col)
+{
+	int i, j, mass=0;
+	for (i=-AVERAGING_RADIUS; i<=AVERAGING_RADIUS; i++) {
+		for (j=-AVERAGING_RADIUS; j<=AVERAGING_RADIUS; j++) {
+			mass+=get_mass(FIELD(i+row, j+col));
+		}
+	}
+	return mass;
+}
+
+void save(int iter)
+{
+	int i, j;
+	char file_name[100];
+	FILE *fl;
+
+	sprintf(file_name, "density/%06d.txt", iter);
+	fl=fopen(file_name, "w");
+
+	for (j=AVERAGING_RADIUS; j<WIDTH-AVERAGING_RADIUS-1; j++) {
+		double density=0;
+		for (i=AVERAGING_RADIUS; i<HEIGHT-AVERAGING_RADIUS-1; i++) {
+			density+=sum_mass(i, j);
+		}
+		density/=HEIGHT*AVERAGING_RADIUS*AVERAGING_RADIUS;
+		fprintf(fl, "%d\t%lf\n", j, density);
+	}
+
+	fclose(fl);
+}
+
+void swap_buffers()
+{
+	char* tmp;
+	tmp=field;
+	field=new_field;
+	new_field=tmp;
+}
+
+void move()
+{
+	int i;
+	for (i=0; i<HEIGHT; i++) {
+		int j;
+		for (j=0; j<WIDTH; j++) {
+			char cell=FIELD(i, j)&0xf0; // keep rest mass
+			
+			if (FIELD(i, (j+WIDTH-1)%WIDTH)&1) cell+=1; // right
+			if (FIELD(i, (j+1)%WIDTH)&4) cell+=4; // left
+			if (FIELD((i+HEIGHT-1)%HEIGHT, j)&8) cell+=8; // up
+			if (FIELD((i+1)%HEIGHT, j)&2) cell+=2; // down
+
+			NEW_FIELD(i, j)=cell;
+		}
+	}
+
+	swap_buffers();
+}
+
+int main()
+{
+	int iter;
+	struct timespec t0, t1;
+
+	field=(char*)malloc(WIDTH*HEIGHT);
+	new_field=(char*)malloc(WIDTH*HEIGHT);
+
+	init();
+
+	clock_gettime(CLOCK_MONOTONIC_RAW, &t0);
+
+	for (iter=0; iter<ITERS; iter++) {
+		int i;
+
+		sources(iter);
+
+		for(i=0; i<HEIGHT; i++) {
+			int j;
+			for (j=0; j<WIDTH; j++) {
+				FIELD(i, j)=collide(FIELD(i, j));
+			}
+		}
+
+		move();
+
+		if (iter%SAVE_PERIOD==0) {
+			save(iter);
+		}
+	}
+
+	clock_gettime(CLOCK_MONOTONIC_RAW, &t1);
+
+	printf("Done %d iters, %lf MCell/sec\n", ITERS,
+			ITERS*HEIGHT*WIDTH/
+			((t1.tv_sec-t0.tv_sec)+0.000000001*(t1.tv_nsec-t0.tv_nsec))
+			/1000000);
+
+	return 0;
+}

image_gen.sh

+#!/bin/sh
+
+for x in `ls $1/*.txt`; do
+	echo 'set term png
+set yrange[0:20]
+set output '\"$x.png\"'
+plot '\"$x\"'' | gnuplot - 2>/dev/null
+done
+
+mencoder mf://$1/*.png -mf fps=10:type=png -ovc lavc -lavcopts vcodec=mpeg4:mbd=2:trell -oac copy -o $1/video.avi