Compare commits

..

5 Commits

Author SHA1 Message Date
c056e12b7f ex0: theory document 2025-08-26 17:52:26 +02:00
6b5adb7f37 ex0: solve :) 2025-08-25 11:52:35 +02:00
4f01b6101b ex0: style format 2025-08-24 11:18:03 +02:00
3a3a378110 ex0: initial 2025-08-23 23:00:35 +02:00
e9e6749bae nixify dependencies 2025-08-23 22:57:31 +02:00
19 changed files with 9 additions and 565 deletions

3
.gitignore vendored
View File

@@ -54,5 +54,4 @@ dkms.conf
*.zip
*.pdf
*.bmp
!main*
*.png

View File

@@ -1,17 +1,13 @@
build: main.c bitmap.c
gcc main.c bitmap.c -o main -O2
build:
gcc main.c bitmap.c -o main -O3
run: build
run:
make build
./main
show: run
show:
make run
feh after.bmp
convert: run
convert:
magick after.bmp after.png # for showing image in pdf
zip: before.bmp main.pdf Makefile main.c bitmap.h bitmap.c
zip handin.zip before.bmp main.pdf Makefile main.c bitmap.h bitmap.c
unzip: handin.zip
unzip handin.zip -d handin

Binary file not shown.

Before

Width:  |  Height:  |  Size: 60 MiB

After

Width:  |  Height:  |  Size: 60 MiB

BIN
exercise0/before.bmp Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 15 MiB

Binary file not shown.

Binary file not shown.

View File

@@ -1,3 +0,0 @@
data
images
*.mp4

View File

@@ -1,19 +0,0 @@
CFLAGS+= -std=c99 -O2 -Wall -Wextra
LDLIBS+= -lm
TARGETS=wave_1d
IMAGES=$(shell find data -type f | sed s/\\.dat/.png/g | sed s/data/images/g )
.PHONY: all clean dirs plot movie
all: dirs ${TARGETS}
dirs:
mkdir -p data images
plot: ${IMAGES}
images/%.png: data/%.dat
./plot_image.sh $<
movie: ${IMAGES}
ffmpeg -y -an -i images/%5d.png -vcodec libx264 -pix_fmt yuv420p -profile:v baseline -level 3 -r 12 wave.mp4
clean:
-rm -fr ${TARGETS} data images wave.mp4
zip: theory.pdf Makefile wave_1d.c
zip handin.zip theory.pdf Makefile wave_1d.c
unzip: handin.zip
unzip handin.zip -d handin

View File

@@ -1,4 +0,0 @@
* make : builds the executable 'wave_1d'
* ./wave\_1d : fills the 'data/' directory with stored time steps
* make plot : converts saved time steps to png files under 'images/', using gnuplot. Runs faster if launched with e.g. 4 threads (make -j4 plot).
* make movie : converts collection of png files under 'images' into an mp4 movie file

View File

@@ -1,10 +0,0 @@
#! /usr/bin/env bash
SIZE=1024
DATAFILE=$1
IMAGEFILE=`echo $1 | sed s/dat$/png/ | sed s/data/images/`
cat <<END_OF_SCRIPT | gnuplot -
set term png
set output "$IMAGEFILE"
set yrange[-1:1]
plot "$DATAFILE" binary array=${SIZE} format='%double' with lines
END_OF_SCRIPT

View File

@@ -1,70 +0,0 @@
#import "@preview/cetz:0.3.2";
#import "@preview/cetz-plot:0.1.1": plot
#import "@preview/physica:0.9.4": *
#import "@preview/plotsy-3d:0.1.0": plot-3d-parametric-surface
#import "@preview/fletcher:0.5.4" as fletcher: diagram, edge, node
#set page(paper: "a4", margin: (x: 2.6cm, y: 2.8cm), numbering: "1 : 1")
#set par(justify: true, leading: 0.52em)
#let FONT_SIZE = 18pt;
#set text(font: "FreeSerif", size: FONT_SIZE, lang: "us")
#show math.equation: set text(font: "Euler Math", size: (FONT_SIZE * 1.0), lang: "en")
#set heading(numbering: none)
#show heading.where(level: 1): it => {
rect(inset: FONT_SIZE / 2)[#it]
}
// Dracula palette
#let bg = rgb("#282a36")
#let fg = rgb("#f8f8f2")
#let sec = rgb("#44475a")
#align(center)[
#text(size: FONT_SIZE * 2, weight: "bold")[#underline[exercise 0]]
]
these are my solutions to exercise set 1 of TDT4200.
this document was created using
#link("https://typst.app/")[#text(blue.darken(5%))[typst]].
#v(42pt)
#outline(title: none)
#v(42pt)
= macros
macros define programmer-friendly syntax that is preprocessed at compile time,
thus incurring no performance penality. they constitute a fundamental part in
meta programming (see lisp). excessive use of macros may obfuscate semantics,
but can often be used to make the code clearer to read. in this case, it helps
simplify the indexing of the buffers for easier to read calculations.
= boundary conditions
the boundary condition could be changed from reflection to wrapping around. this
can be achieved by setting the left boundary to be the right-most point, and the
right boundary to be the left-most point.
= no memory
if you don't allocate memory in T1, the buffers will be unallocated and you will
get indexing errors as you try to access the buffers using the predefined
macros.
basically, you segfault.
= `float const *a` vs `float *const b`
- `float const *a` is a constant pointer to some memory address, assigned to the
constant name `a`. since `a` is constant, it cannot be reassigned.
- `float *const b` is a pointer to some constant value, assigned to the
variable `b`. since `b` points to a constant, the dereferenced value cannot be
changed.
these are not the same.

Binary file not shown.

View File

@@ -1,120 +0,0 @@
#define _XOPEN_SOURCE 600
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// Option to change numerical precision.
typedef int64_t int_t;
typedef double real_t;
// Simulation parameters: size, step count, and how often to save the state.
const int_t
N = 1024,
max_iteration = 4000,
snapshot_freq = 10;
// Wave equation parameters, time step is derived from the space step.
const real_t
c = 1.0,
dx = 1.0;
real_t
dt;
// Buffers for three time steps, indexed with 2 ghost points for the boundary.
real_t
*buffers[3] = { NULL, NULL, NULL };
#define U_prv(i) buffers[0][(i) + 1]
#define U(i) buffers[1][(i) + 1]
#define U_nxt(i) buffers[2][(i) + 1]
// Save the present time step in a numbered file under 'data/'.
void domain_save(int_t step) {
char filename[256];
sprintf(filename, "data/%.5ld.dat", step);
FILE *out = fopen(filename, "wb");
fwrite(&U(0), sizeof(real_t), N, out);
fclose(out);
}
// TASK: T1
// Set up our three buffers, fill two with an initial cosine wave,
// and set the time step.
void domain_initialize(void) {
// allocate memory
for (int i = 0; i < 3; i++)
buffers[i] = malloc((N + 2) * sizeof(real_t));
// apply initial condition
real_t x = 0.0f;
for (int i = 0; i < N; i++) {
U_prv(i) = cos(2 * M_PI * x);
U(i) = U_prv(i);
U_nxt(i) = 0.0f;
x += (real_t)dx / (real_t)N;
}
// set the time-step dt according to CFL
dt = dx / c;
}
// TASK T2:
// Return the memory to the OS.
void domain_finalize(void) {
for (int i = 0; i < 3; i++) {
free(buffers[i]);
buffers[i] = NULL;
}
}
// TASK: T3
// Rotate the time step buffers.
void rotate_buffers(void) {
real_t *temp = buffers[0];
buffers[0] = buffers[1];
buffers[1] = buffers[2];
buffers[2] = temp;
}
// TASK: T4
// Derive step t+1 from steps t and t-1.
void time_step(void) {
for (int i = 0; i < N; i++) {
U_nxt(i) = -U_prv(i) + 2 * U(i) + dt * dt * c * c / (dx * dx) * (U(i - 1) + U(i + 1) - 2 * U(i));
}
}
// TASK: T5
// Neumann (reflective) boundary condition.
void boundary(void) {
U_prv(-1) = U_prv(1);
U(-1) = U(1);
U_nxt(-1) = U_nxt(1);
U_prv(N) = U_prv(N - 2);
U(N) = U(N - 2);
U_nxt(N) = U_nxt(N - 2);
}
// TASK: T6
// Main time integration.
void simulate(void) {
int_t iteration = 0;
while (iteration < max_iteration) {
boundary();
time_step();
rotate_buffers();
domain_save(iteration / snapshot_freq);
iteration++;
}
}
int main(void) {
domain_initialize();
simulate();
domain_finalize();
exit(EXIT_SUCCESS);
}

View File

@@ -1 +0,0 @@
out/

View File

@@ -1,33 +0,0 @@
CC := gcc
CFLAGS := -Wall -Wextra -std=c17
SRC := mandel_mpi.c
TARGET := $(SRC:.c=)
OUTDIR := out
OUT := $(OUTDIR)/$(TARGET)
ARGS := 1
NBENCH := 1
NPROC := 3
.PHONY: all clean run parallel
all: clean parallel show
$(OUT): $(SRC)
$(CC) $(CFLAGS) -o $(OUT) $<
run: $(OUT)
cd $(OUTDIR) && ./$(TARGET) $(ARGS)
clean:
rm -rf $(OUTDIR)/*
time: $(OUT)
python3 bench.py './$(OUT) 0' $(NBENCH)
parallel: $(SRC)
mpicc -o $(OUT) $(SRC)
cd $(OUTDIR) && mpirun -np $(NPROC) $(TARGET) $(ARGS)
show: $(OUTDIR)/mandel.bmp
feh $<

View File

@@ -1,24 +0,0 @@
#!/usr/bin/env python3
from statistics import mean
from subprocess import DEVNULL, run
from sys import argv, exit
from time import perf_counter
def benchmark(cmd, N=10):
times = []
for _ in range(N):
start = perf_counter()
run(cmd, shell=True, stdout=DEVNULL, stderr=DEVNULL)
times.append(perf_counter() - start)
return mean(times)
if __name__ == "__main__":
if len(argv) < 2:
print("usage: python3 bench.py 'cmd' [N]")
exit(1)
cmd = argv[1]
N = int(argv[2]) if len(argv) > 2 else 10
print(f"{benchmark(cmd, N):.4f}")

View File

@@ -1,100 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define XSIZE 2560
#define YSIZE 2048
#define MAXITER 255
double xleft=-2.01;
double xright=1;
double yupper,ylower;
double ycenter=1e-6;
double step;
int pixel[XSIZE*YSIZE];
#define PIXEL(i,j) ((i)+(j)*XSIZE)
typedef struct {
double real,imag;
} complex_t;
void calculate() {
for(int i=0;i<XSIZE;i++) {
for(int j=0;j<YSIZE;j++) {
/* Calculate the number of iterations until divergence for each pixel.
If divergence never happens, return MAXITER */
complex_t c,z,temp;
int iter=0;
c.real = (xleft + step*i);
c.imag = (ylower + step*j);
z = c;
while(z.real*z.real + z.imag*z.imag < 4) {
temp.real = z.real*z.real - z.imag*z.imag + c.real;
temp.imag = 2*z.real*z.imag + c.imag;
z = temp;
if(++iter==MAXITER) break;
}
pixel[PIXEL(i,j)]=iter;
}
}
}
typedef unsigned char uchar;
/* save 24-bits bmp file, buffer must be in bmp format: upside-down */
void savebmp(char *name,uchar *buffer,int x,int y) {
FILE *f=fopen(name,"wb");
if(!f) {
printf("Error writing image to disk.\n");
return;
}
unsigned int size=x*y*3+54;
uchar header[54]={'B','M',size&255,(size>>8)&255,(size>>16)&255,size>>24,0,
0,0,0,54,0,0,0,40,0,0,0,x&255,x>>8,0,0,y&255,y>>8,0,0,1,0,24,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
fwrite(header,1,54,f);
fwrite(buffer,1,XSIZE*YSIZE*3,f);
fclose(f);
}
/* given iteration number, set a colour */
void fancycolour(uchar *p,int iter) {
if(iter==MAXITER);
else if(iter<8) { p[0]=128+iter*16; p[1]=p[2]=0; }
else if(iter<24) { p[0]=255; p[1]=p[2]=(iter-8)*16; }
else if(iter<160) { p[0]=p[1]=255-(iter-24)*2; p[2]=255; }
else { p[0]=p[1]=(iter-160)*2; p[2]=255-(iter-160)*2; }
}
int main(int argc,char **argv) {
if(argc==1) {
puts("Usage: MANDEL n");
puts("n decides whether image should be written to disk (1=yes, 0=no)");
return 0;
}
/* Calculate the range in the y-axis such that we preserve the
aspect ratio */
step=(xright-xleft)/XSIZE;
yupper=ycenter+(step*YSIZE)/2;
ylower=ycenter-(step*YSIZE)/2;
calculate();
if(strtol(argv[1],NULL,10)!=0) {
/* create nice image from iteration counts. take care to create it upside
down (bmp format) */
unsigned char *buffer=calloc(XSIZE*YSIZE*3,1);
for(int i=0;i<XSIZE;i++) {
for(int j=0;j<YSIZE;j++) {
int p=((YSIZE-j-1)*XSIZE+i)*3;
fancycolour(buffer+p,pixel[PIXEL(i,j)]);
}
}
/* write image to disk */
savebmp("mandel2.bmp",buffer,XSIZE,YSIZE);
}
return 0;
}

View File

@@ -1,158 +0,0 @@
#include <math.h>
#include <mpi.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define XSIZE 2560
#define YSIZE 2048
#define MAXITER 255
double xleft = -2.01;
double xright = 1;
double yupper, ylower;
double ycenter = 1e-6;
double step;
int pixel[XSIZE * YSIZE];
#define PIXEL(i, j) ((i) + (j) * XSIZE)
typedef struct {
double real, imag;
} complex_t;
int rank, size;
bool save = false;
void calculate() {
for (int i = 0; i < XSIZE; i++) {
for (int j = 0; j < YSIZE; j++) {
/* Calculate the number of iterations until divergence for each pixel.
If divergence never happens, return MAXITER */
complex_t c, z, temp;
int iter = 0;
c.real = (xleft + step * i);
c.imag = (ylower + step * j);
z = c;
while (z.real * z.real + z.imag * z.imag < 4) {
temp.real = z.real * z.real - z.imag * z.imag + c.real;
temp.imag = 2 * z.real * z.imag + c.imag;
z = temp;
if (++iter == MAXITER)
break;
}
pixel[PIXEL(i, j)] = iter;
}
}
}
typedef unsigned char uchar;
/* save 24-bits bmp file, buffer must be in bmp format: upside-down */
void savebmp(char *name, uchar *buffer, int x, int y) {
FILE *f = fopen(name, "wb");
if (!f) {
printf("Error writing image to disk.\n");
return;
}
unsigned int size = x * y * 3 + 54;
uchar header[54] = { 'B', 'M', size & 255, (size >> 8) & 255, (size >> 16) & 255, size >> 24, 0,
0, 0, 0, 54, 0, 0, 0, 40, 0, 0, 0, x & 255, x >> 8, 0, 0, y & 255, y >> 8, 0, 0, 1, 0, 24, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
fwrite(header, 1, 54, f);
fwrite(buffer, 1, XSIZE * YSIZE * 3, f);
fclose(f);
}
/* given iteration number, set a colour */
void fancycolour(uchar *p, int iter) {
if (iter == MAXITER)
;
else if (iter < 8) {
p[0] = 128 + iter * 16;
p[1] = p[2] = 0;
} else if (iter < 24) {
p[0] = 255;
p[1] = p[2] = (iter - 8) * 16;
} else if (iter < 160) {
p[0] = p[1] = 255 - (iter - 24) * 2;
p[2] = 255;
} else {
p[0] = p[1] = (iter - 160) * 2;
p[2] = 255 - (iter - 160) * 2;
}
}
int main(int argc, char **argv) {
if (argc <= 1) {
printf("usage: mandel_mpi [y/n]\n");
printf("y/n denotes if you want to save the resulting\n");
printf("calculations to mandel.bmp or not.\n");
return 1;
} else {
char arg = argv[1][0];
save = arg == 'y' || arg == '1';
}
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
#undef YSIZE
#define YSIZE (2048 / size)
step = (xright - xleft) / XSIZE;
double old_yupper = ycenter + (step * YSIZE * size) / 2;
double old_ylower = ycenter - (step * YSIZE * size) / 2;
double delta = (old_yupper - old_ylower) / size;
ylower = old_ylower + rank * delta;
yupper = old_ylower + (rank + 1) * delta;
ycenter = (ylower + yupper) / 2;
calculate();
printf("after calculate %d\n", rank);
const int bufsize = XSIZE * YSIZE * 3;
unsigned char *buffer = calloc(bufsize, 1);
for (int i = 0; i < XSIZE; i++) {
for (int j = 0; j < YSIZE; j++) {
int p = ((YSIZE - j - 1) * XSIZE + i) * 3;
fancycolour(buffer + p, pixel[PIXEL(i, j)]);
}
}
printf("after buffer %d\n", rank);
if (rank != 0)
MPI_Send(buffer, bufsize, MPI_UNSIGNED_CHAR, 0, 0, MPI_COMM_WORLD);
printf("hello again from %d\n", rank);
// use process 0 as central process: gather calulations and output image
if (rank == 0) {
unsigned char *total = calloc(bufsize, size);
memcpy(total + (size - 1) * bufsize, buffer, bufsize);
printf("after memcpy 0\n");
for (int i = 1; i < size; i++) {
MPI_Recv(total + (size - 1 - i) * bufsize, bufsize, // write reverse order
MPI_UNSIGNED_CHAR, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
printf("after direct receive %d\n", i);
}
printf("after loop\n");
#undef YSIZE
#define YSIZE 2048
if (save)
savebmp("mandel.bmp", total, XSIZE, YSIZE);
printf("after save %d\n", save);
}
MPI_Finalize();
return 0;
}

View File

@@ -16,17 +16,8 @@
buildInputs = with pkgs; [
gcc
gnumake
imagemagick_light
typst
typstyle
zip
unzip
feh
# exercise 1: 1D wave equation plotting utils
gnuplot
ffmpeg
mpi
python3
feh
];
shellHook = ''
echo welcome!