Added code for hands-on

_quarto.yml

@@ -52,7 +52,7 @@ website:
           - "exercises/user-experience.qmd"
           # - "exercises/testing.qmd"
           # - "exercises/git2.qmd"
-          - "exercises/parallelism.qmd"
+          - "exercises/parallelism/parallelism.qmd"
           # - "exercises/hardware.qmd"
           # - "exercises/file_and_data_systems.qmd"
           # - "exercises/memory_hierarchies.qmd"

exercises/parallelism/main.c

+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+
+#ifdef _OPENMP
+#include <omp.h>
+#else
+#ifdef WRITE_DECOMP
+#error You cannot write the decomposition without OpenMP.
+#endif
+#endif
+
+//compile using (for example):
+//gcc -lm -fopenmp -g main.c
+
+//definition of constants
+//some values are given in meter
+//some in seconds
+//
+#define x_0 10000
+#define delta_x 1
+#define coeff_sech 0.01
+//depth of sea at deepest point
+#define b_0 1000
+//gravitational acceleration
+#define g 10
+#define delta_t 0.001
+
+
+double *h;
+double *h1;
+double *h2;
+double *u;
+double *u1;
+double *u2;
+
+int *point_to_thread;
+
+long count_points () {
+	return (int)(2*x_0/delta_x);
+}
+
+double get_x(long i){
+	return i*delta_x - x_0;
+}
+
+double waterdepth(long i){
+	double x = get_x(i);
+	return b_0*(x*x/(x_0*x_0)-1);
+}
+
+void initialise_h(double *h){
+	double x, tmp;
+	long i;
+	for(i=0; i < count_points(); i++){
+		x=get_x(i);
+		tmp = 2/(exp(coeff_sech*x)+exp(-coeff_sech*x));
+		if (tmp > 1e-10)
+			h[i] = tmp*tmp;
+		else
+			h[i] = 0;
+	}
+}
+
+void initialise_u(double *u){
+	long i;
+	for(i=0;i<count_points();i++){
+		u[i]=0;
+	}
+}
+
+
+
+void print_values(FILE *data, double *h, double *u){
+	long i;
+	fprintf(data, "#x\th\tu\tb\n");
+	for(i=0;i<count_points();i++){
+		fprintf(data, "%g\t%g\t%g\t%g\n", get_x(i), h[i], u[i], waterdepth(i));
+	}
+}
+
+double lu(double *h, double *u, long i){
+	return u[i]*(u[i+1]-u[i-1])+g*(h[i+1]-h[i-1]);
+}
+
+double lh(double *h, double *u, long i){
+	return u[i]*(h[i+1]-h[i-1])-u[i]*(waterdepth(i+1)-waterdepth(i-1))+(h[i]-waterdepth(i))*(u[i+1]-u[i-1]);
+}
+
+
+//sets h1,u1 using values (startconditions) from u2/h2
+void euler(double *h1, double *h2, double *u1, double *u2){
+	long i;
+	long npoints = count_points();
+	u1[0]=0;		//boundary value
+	h1[0]=0;		//boundary value
+	u1[npoints-1]=0;
+	h1[npoints-1]=0;
+	for(i=1;i<npoints-1;i++){
+		u1[i]=u2[i]-(delta_t/(2.*delta_x))*lu(h2,u2,i);
+		h1[i]=h2[i]-(delta_t/(2.*delta_x))*lh(h2,u2,i);
+	}
+}
+
+long findmax(double *h, long oldmaxi,long outputmth){
+	long i, maxi=oldmaxi-10;
+	for(i=oldmaxi-10;(i< oldmaxi+(int)(outputmth/delta_x)+30) && (i<count_points()) ;i++){
+		if(h[i]>h[maxi]){
+			maxi=i;
+		}
+	}
+	return maxi;
+}
+
+double maxval(double *h) {
+	long npoints = count_points();
+	long i;
+	double hmax = -10;
+// HOMEWORK: Parallelise this loop!
+	for(i = 0; i < npoints; i++){
+		if (h[i] > hmax)
+			hmax = h[i];
+	}
+	return hmax;
+}
+
+
+double ng_derivative(double xi, double xi_1, double deltax){
+	return (xi-xi_1)/deltax;
+}
+
+
+void siginthandler(){
+	print_values(stdout,h,u);
+	exit(0);
+}
+
+#ifdef _OPENMP
+void save_decomp(long i) {
+	point_to_thread[i] = omp_get_thread_num();
+}
+
+void print_decomp() {
+	int i,j;
+	int npoints = count_points();
+	for (i=0; i<omp_get_num_threads(); i++) {
+		#pragma omp barrier
+		if (i != omp_get_thread_num())
+			continue;
+		printf("%d\t", i);
+		int first=-1;
+		for (j=0; j<npoints+1; j++) {
+			if (j < npoints && point_to_thread[j] == i) {
+				if (first == -1)
+					first = j;
+			} else {
+				if (first == -1)
+					continue;
+				if (first == j-1)
+					printf("%d ", j);
+				else {
+					printf("%d-%d ", first, j-1);
+				}
+				first = -1;
+			}
+		}
+		printf("\n");
+	}
+}
+#endif
+
+
+void swap(double **var, double **var1, double **var2) {
+	double *temp;
+
+	temp  = *var2;
+	*var2 = *var1;
+	*var1 = *var;
+	*var  = temp;
+}
+
+void set_boundaries(double *var) {
+	long npoints = count_points();
+	var[0]         = 0;
+	var[npoints-1] = 0;
+}
+
+int main(int argc, char* argv){
+	long i;
+	long npoints = count_points();
+	h  = malloc(sizeof(double)*npoints);
+	h1 = malloc(sizeof(double)*npoints);
+	h2 = malloc(sizeof(double)*npoints);
+
+	u  = malloc(sizeof(double)*npoints);
+	u1 = malloc(sizeof(double)*npoints);
+	u2 = malloc(sizeof(double)*npoints);
+#ifdef WRITE_DECOMP
+	point_to_thread = malloc(sizeof(int)*count_points());
+#endif
+
+	if(h==NULL || h1==NULL || h2==NULL || u==NULL || u1==NULL || u2==NULL){
+		fprintf(stderr, "Insufficient memory");
+		exit(23);
+	}
+#if 0
+	signal(SIGINT, siginthandler);
+#endif
+
+	long maxi=((npoints)/2)+5;
+	long oldmaxi;
+	long oldmaxi_output=(npoints)/2;
+	int maximth=20;
+	double t=0;
+	long n=0;
+	int outputmth=1000; //write every mth set of values to stderr
+	int dataqth=200; //write every qth set of values to output#
+	char file[60]; //char pointer for filename
+	FILE* out;
+	double *temp;
+	int first = 1;
+
+	initialise_u(u2);
+	initialise_h(h2);
+	euler(h1,h2,u1,u2);
+#if 0
+        out = fopen("t-1.dat","w");
+        print_values(out,h1,u1);
+        out = fopen("t-2.dat","w");
+        print_values(out,h2,u2);
+#endif
+	time_t t_before = time(NULL);
+	fprintf(stderr, "#time	h_max\n");
+
+	//while ( t < 150 ){
+	while ( t < 50 ){
+		if (!first) {
+			swap(&h, &h1, &h2);
+			swap(&u, &u1, &u2);
+		} else
+			first = 0;
+
+		set_boundaries(u);
+		set_boundaries(h);
+
+		#pragma omp parallel for // schedule(dynamic, 10)
+		for(i=0;i<npoints;i++){
+#ifdef WRITE_DECOMP
+			save_decomp(i);
+#endif
+			if (i == 0 || i == npoints-1)
+				continue;
+			u[i]=u2[i]-(delta_t/delta_x)*lu(h1,u1,i);
+			h[i]=h2[i]-(delta_t/delta_x)*lh(h1,u1,i);
+		}
+
+		t+=delta_t;
+		if((n%maximth)==0){
+
+			oldmaxi=maxi;
+			maxi=findmax(h,oldmaxi,maximth);
+			if((n%outputmth)==0){
+				double maxh = maxval(h);
+				fprintf(stderr, "%g\t%g\n", t, maxh);
+
+#if 0
+				fprintf(stderr, "%g \t",t);
+				fprintf(stderr, "%g\t%g \t%g\n",get_x(maxi),h[maxi],ng_derivative(get_x(maxi),get_x(oldmaxi_output),outputmth*delta_t));
+				oldmaxi_output=maxi;
+#endif
+			}
+		}
+
+
+#if 0
+                  if((n%dataqth)==0){
+                        sprintf(file,"output%04ld.dat",n/dataqth);
+                        out = fopen(file,"w");
+                        print_values(out,h,u);
+                  }
+#endif
+		n++;
+	}
+	time_t t_after = time(NULL);
+	double duration = difftime(t_after, t_before);
+	printf("Run took %g seconds.\n", duration);
+
+	//print_values(stdout,h,u);
+#ifdef WRITE_DECOMP
+#pragma omp parallel
+	print_decomp();
+#endif
+
+	return 0;
+}
+