Commit 8ba76279 authored by Vladislav Perepelkin's avatar Vladislav Perepelkin
Browse files

improved image generation

parent 8023d0a8
#include <assert.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
......@@ -15,6 +16,7 @@
#define RP4 128
char *field, *new_field;
int ensemble;
void fill(int i0, int j0, int i1, int j1,
double r, double d, double l, double u,
......@@ -75,19 +77,65 @@ void sum_mass(int row, int col, int* rest, int* move)
}
}
double* load(char const* path)
{
int j;
FILE *fl;
double *data=(double*)malloc(
(WIDTH-2*AVERAGING_RADIUS-1)*sizeof(double)
);
fl=fopen(path, "r");
for(j=AVERAGING_RADIUS; j<WIDTH-AVERAGING_RADIUS-1; j++) {
int row;
double val;
if (fscanf(fl, "%d\t%lf\n", &row, &val)!=2 || row!=j) {
fprintf(stderr, "ERROR: File '%s' has illegal format",
path);
abort();
} else {
data[j-AVERAGING_RADIUS]=val;
}
}
fclose(fl);
return data;
}
void write(FILE *fl, int j, int ensemble, double *old, double val)
{
assert(ensemble? !!old: !old);
fprintf(fl, "%d\t%lf\n", j,
ensemble?
(val+ensemble*old[j-AVERAGING_RADIUS])/(ensemble+1)
: val
);
}
void save(int iter)
{
int i, j;
char file_name[100];
char density_path[100], move_path[100], rest_path[100];
FILE *fl_density, *fl_move, *fl_rest;
double square=(AVERAGING_RADIUS*2+1)*(AVERAGING_RADIUS*2+1);
double *old_density=0, *old_move=0, *old_rest=0;
sprintf(density_path, "density/%06d.txt", iter);
sprintf(move_path, "move/%06d.txt", iter);
sprintf(rest_path, "rest/%06d.txt", iter);
if (ensemble>0) {
old_density=load(density_path);
old_move=load(move_path);
old_rest=load(rest_path);
}
sprintf(file_name, "density/%06d.txt", iter);
fl_density=fopen(file_name, "w");
sprintf(file_name, "move/%06d.txt", iter);
fl_move=fopen(file_name, "w");
sprintf(file_name, "rest/%06d.txt", iter);
fl_rest=fopen(file_name, "w");
fl_density=fopen(density_path, "w");
fl_move=fopen(move_path, "w");
fl_rest=fopen(rest_path, "w");
for (j=AVERAGING_RADIUS; j<WIDTH-AVERAGING_RADIUS-1; j++) {
int rest=0, move=0;
......@@ -95,14 +143,22 @@ void save(int iter)
sum_mass(i, j, &rest, &move);
}
fprintf(fl_density, "%d\t%lf\n", j, 1.0*(rest+move)/(HEIGHT-AVERAGING_RADIUS*2-1)/square);
fprintf(fl_move, "%d\t%lf\n", j, 1.0*move/(HEIGHT-AVERAGING_RADIUS*2-1)/square);
fprintf(fl_rest, "%d\t%lf\n", j, 1.0*rest/(HEIGHT-AVERAGING_RADIUS*2-1)/square);
write(fl_density, j, ensemble, old_density,
1.0*(rest+move)/(HEIGHT-AVERAGING_RADIUS*2-1)/square);
write(fl_move, j, ensemble, old_move,
1.0*move/(HEIGHT-AVERAGING_RADIUS*2-1)/square);
write(fl_rest, j, ensemble, old_rest,
1.0*rest/(HEIGHT-AVERAGING_RADIUS*2-1)/square);
}
fclose(fl_density);
fclose(fl_rest);
fclose(fl_move);
free(old_density);
free(old_move);
free(old_rest);
}
void swap_buffers()
......@@ -120,8 +176,6 @@ int main()
field=(char*)malloc(WIDTH*HEIGHT);
new_field=(char*)malloc(WIDTH*HEIGHT);
init();
clock_gettime(CLOCK_MONOTONIC_RAW, &t0);
#pragma omp parallel num_threads(THREADS)
......@@ -130,49 +184,62 @@ int main()
#pragma omp single
printf("Running %d threads\n", omp_get_num_threads());
for (iter=0; iter<ITERS; iter++) {
int i;
for (ensemble=0; ensemble<ENSEMBLE_SIZE;)
{
#pragma omp single
sources(iter);
printf("ENSEMBLE: %d\n", ensemble);
init();
for (iter=0; iter<ITERS; iter++) {
int i;
#pragma omp for private(i) schedule(static)
for(i=0; i<HEIGHT; i++) {
int j;
for (j=0; j<WIDTH; j++) {
FIELD(i, j)=collide(FIELD(i, j));
#pragma omp single
sources(iter);
#pragma omp for private(i) schedule(static)
for(i=0; i<HEIGHT; i++) {
int j;
for (j=0; j<WIDTH; j++) {
FIELD(i, j)=collide(FIELD(i, j));
}
}
}
#pragma omp barrier
#pragma omp barrier
#pragma omp for private(i) schedule(static)
for (i=0; i<HEIGHT; i++) {
int j;
for (j=0; j<WIDTH; j++) {
char cell=FIELD(i, j)&0xf0; // keep rest mass
// propogate particles
if (FIELD(i, (j+WIDTH-1)%WIDTH)&RIGHT) cell+=RIGHT;
if (FIELD(i, (j+1)%WIDTH)&LEFT) cell+=LEFT;
if (FIELD((i+HEIGHT-1)%HEIGHT, j)&UP) cell+=UP;
if (FIELD((i+1)%HEIGHT, j)&DOWN) cell+=DOWN;
NEW_FIELD(i, j)=cell;
}
}
#pragma omp for private(i) schedule(static)
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)&RIGHT) cell+=RIGHT; // right
if (FIELD(i, (j+1)%WIDTH)&LEFT) cell+=LEFT; // left
if (FIELD((i+HEIGHT-1)%HEIGHT, j)&UP) cell+=UP; // up
if (FIELD((i+1)%HEIGHT, j)&DOWN) cell+=DOWN; // down
#pragma omp barrier
#pragma omp single
{
swap_buffers();
NEW_FIELD(i, j)=cell;
if (iter%SAVE_PERIOD==0) {
save(iter);
printf("%d/%d\r", iter, ITERS);
fflush(stdout);
}
}
#pragma omp barrier
}
#pragma omp barrier
#pragma omp single
{
swap_buffers();
if (iter%SAVE_PERIOD==0) {
save(iter);
printf("%d/%d\r", iter, ITERS);
fflush(stdout);
}
}
#pragma omp barrier
ensemble++;
}
}
......
#!/bin/sh
zmin=3
zmax=5
zmin=3.2
zmax=4.7
for x in `ls $1/*.txt`; do
echo 'set term png
set yrange['$zmin':'$zmax']
set output '\"$x.png\"'
plot '\"$x\"'' | gnuplot - 2>/dev/null
set grid ytics lw 1 lt 0
set grid xtics lw 1 lt 0
unset key
plot '\"$x\"' with lines lw 4' | 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 $2 >/dev/null 2>/dev/null
tmp_script=`tempfile`
echo 'set term png
set yrange['$zmin':'$zmax']
set output '\"$x.png\"'
set grid ytics lw 1 lt 0
set grid xtics lw 1 lt 0
unset key' > $tmp_script
first=true
for x in `ls $1/*.txt`; do
$first && echo "plot '$x' with lines lw 4\\" >> $tmp_script
$first || echo ", '$x' with lines lw 4\\" >> $tmp_script
first=false
done
gnuplot $tmp_script 2>/dev/null
rm $tmp_script
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment