tFirst commit - werner - cellular automata simulation of wind-driven sand transport
(HTM) git clone git://src.adamsgaard.dk/werner
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit d238e187ed68520db7a7d73e1fb9b011ec1d8c53
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Sun, 21 Apr 2013 11:08:31 +0200
First commit
Diffstat:
A Makefile | 35 +++++++++++++++++++++++++++++++
A README.rst | 31 +++++++++++++++++++++++++++++++
A initrnd.c | 29 +++++++++++++++++++++++++++++
A iterate.c | 35 +++++++++++++++++++++++++++++++
A out/plotmatrix.gp | 20 ++++++++++++++++++++
A printmat.c | 26 ++++++++++++++++++++++++++
A werner.c | 318 +++++++++++++++++++++++++++++++
A werner.h | 23 +++++++++++++++++++++++
A wernerparams.h | 11 +++++++++++
9 files changed, 528 insertions(+), 0 deletions(-)
---
(DIR) diff --git a/Makefile b/Makefile
t@@ -0,0 +1,35 @@
+CFLAGS=-g -Wall -O3 -std=gnu11 -pg
+LDLIBS=`pkg-config --libs gsl`
+LDFLAGS=-pg
+
+default: plots
+
+initrnd: initrnd.o werner.o *.h
+
+iterate: iterate.o werner.o *.h
+
+printmat: printmat.o werner.o *.h
+
+display: matrix.png
+ display $<
+
+matrixes: initrnd iterate printmat Makefile
+ ./initrnd > tmp.mat
+ for i in {001..400}; do cat tmp.mat | ./iterate 1 > tmp_new.mat && cat tmp_new.mat | ./printmat > out/matrix.$$i.txt && mv tmp_new.mat tmp.mat; done
+
+plots: matrixes
+ cd out && for f in *.txt; do gnuplot -e "matrixfile='$$f'" plotmatrix.gp; done
+
+view: plots
+ feh out/*.png
+
+profile-iter: initrnd iterate
+ ./initrnd > tmp.mat
+ ./iterate 10 > tmp_new.mat
+
+clean:
+ $(RM) initrnd iterate printmat
+ $(RM) *.o
+ $(RM) tmp.mat
+ $(RM) out/*.txt
+ $(RM) out/*.png
(DIR) diff --git a/README.rst b/README.rst
t@@ -0,0 +1,31 @@
+Werner
+======
+
+Discrete simulation of wind driven (
+`aeolian <https://en.wikipedia.org/wiki/Aeolian_processes>`_) erosion, transport,
+deposition and avalanching of sand, based on `B.T. Werner's algorithm
+<http://geology.geoscienceworld.org/content/23/12/1107>`_.
+
+The software is licensed under the GNU General Public License, v.3. See
+`license.txt` for more information. These programs are created by Anders
+Damsgaard, `<anders.damsgaard@geo.au.dk>`_.
+
+Requirements
+------------
+The requirements are:
+ * A C compiler toolkit, e.g. the `GNU Compiler Collection
+ <http://gcc.gnu.org/>`_ (GCC)
+ * `GNU Make <https://www.gnu.org/software/make/>`_
+ * `Gnuplot <http://www.gnuplot.info/>`_
+ * `GSL - GNU Scientific Library <http://www.gnuplot.info/>`_
+
+Build and run instructions
+--------------------------
+From the root directory, execute the following command::
+ make
+
+This will build the program, and run an example simulation. Output images are
+placed in the `out/` folder.
+Simulation parameters are changed in `wernerparams.h`. After changing this file,
+run ``make`` again.
+
(DIR) diff --git a/initrnd.c b/initrnd.c
t@@ -0,0 +1,29 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+// see https://www.gnu.org/software/gsl/manual/html_node/Matrices.html
+#include <gsl/gsl_matrix.h>
+
+#include "wernerparams.h"
+#include "werner.h"
+
+
+int main(int argc, char** argv)
+{
+
+ // Allocate matrix Z: number of sand slabs
+ gsl_matrix* Z = gsl_matrix_alloc(rows, cols);
+
+ // Initialize with random values
+ init_rnd_matrix(Z, 1.0, 10.0);
+
+ // Add 100 sand slabs to each cell
+ //gsl_matrix_add_constant(Z, 100.0);
+
+ // Write matrix to stdout
+ gsl_matrix_fprintf(stdout, Z, "%f");
+
+ // End program
+ gsl_matrix_free(Z);
+ return 0;
+}
(DIR) diff --git a/iterate.c b/iterate.c
t@@ -0,0 +1,35 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+// see https://www.gnu.org/software/gsl/manual/html_node/Matrices.html
+#include <gsl/gsl_matrix.h>
+
+#include "wernerparams.h"
+#include "werner.h"
+
+
+int main(int argc, char** argv)
+{
+ // check if the user specified the number of times to iterate
+ int N_c = 1; // Number of cycles, default value
+ if (argc > 1)
+ N_c = atoi(argv[1]);
+
+ // Allocate matrix Z: number of sand slabs
+ gsl_matrix* Z = gsl_matrix_alloc(rows, cols);
+
+ // Read matrix from stdin
+ FILE* fp = fopen("tmp.mat", "r");
+ gsl_matrix_fscanf(fp, Z);
+ //gsl_matrix_fscanf(stdin, Z);
+
+ // Run N_c Werner iterations
+ werner_loop(Z, N_c, d_l, p_ns, p_s);
+
+ // Print result to stdout
+ gsl_matrix_fprintf(stdout, Z, "%f");
+
+ // End program
+ gsl_matrix_free(Z);
+ return 0;
+}
(DIR) diff --git a/out/plotmatrix.gp b/out/plotmatrix.gp
t@@ -0,0 +1,20 @@
+#!/usr/env gnuplot
+
+# plot matrix file. Invoke with:
+# $ gnuplot -e "matrixfile='file.txt'" plotmatrix.gp
+
+set terminal pngcairo
+set output matrixfile.".png"
+
+set size ratio -1
+
+set title "Werner model for aeolian transport and deposition"
+set xlabel "x"
+set ylabel "y"
+set cblabel "sand slabs"
+
+# With interpolation
+set pm3d map interpolate 0,0
+set palette gray
+splot matrixfile matrix
+
(DIR) diff --git a/printmat.c b/printmat.c
t@@ -0,0 +1,26 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+// see https://www.gnu.org/software/gsl/manual/html_node/Matrices.html
+#include <gsl/gsl_matrix.h>
+
+#include "wernerparams.h"
+#include "werner.h"
+
+
+int main(int argc, char** argv)
+{
+
+ // Allocate matrix Z: number of sand slabs
+ gsl_matrix* Z = gsl_matrix_alloc(rows, cols);
+
+ // Read matrix from stdin
+ gsl_matrix_fscanf(stdin, Z);
+
+ // Print formatted result to stdout
+ print_matrix(Z);
+
+ // End program
+ gsl_matrix_free(Z);
+ return 0;
+}
(DIR) diff --git a/werner.c b/werner.c
t@@ -0,0 +1,318 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <gsl/gsl_matrix.h>
+
+// Return a random number in [0;1[
+double rnd(void)
+{
+ return (double)rand()/RAND_MAX;
+}
+
+// Initialize matrix M with random ints from low to high
+void init_rnd_matrix(gsl_matrix* M, double low, double high)
+{
+ int row, col;
+ for (row = 0; row < M->size1; row++)
+ for (col = 0; col < M->size2; col++)
+ gsl_matrix_set(M, row, col,
+ (int)(rnd() * (high-low) + low));
+}
+
+// Write matrix to stdout
+void print_matrix(gsl_matrix* M)
+{
+ int row, col;
+ double val;
+ for (row = 0; row < M->size1; row++) {
+ for (col = 0; col < M->size2; col++) {
+ val = gsl_matrix_get(M, row, col);
+ printf("%f\t", val);
+ }
+ puts(""); // newline
+ }
+}
+
+// Add a slab of sand from the target cell
+void deposit(gsl_matrix* Z, int row, int col)
+{
+ gsl_matrix_set(Z, row, col, gsl_matrix_get(Z, row, col) + 1.0);
+}
+
+// Remove a slab of sand from the target cell
+void erode(gsl_matrix* Z, int row, int col)
+{
+ gsl_matrix_set(Z, row, col, gsl_matrix_get(Z, row, col) - 1.0);
+}
+
+// Find out if cell is in shade (1) or not (0)
+int inshade(
+ gsl_matrix* Z, // matrix with sand slabs
+ int row, // row idx of interest
+ int col) // col idx of interest
+{
+ int js;
+ int i = 1;
+ int shade = 0;
+ int nx = Z->size2;
+
+ while ((i <= nx/4) && (shade == 0)) {
+ js = col - i;
+ if (js < 0) // periodic boundary
+ js += nx;
+ if (gsl_matrix_get(Z, row, js) >= gsl_matrix_get(Z, row, col) + i)
+ shade = 1;
+ i++;
+ }
+ return shade;
+}
+
+void find_max_slope_neighbor_ero(
+ gsl_matrix* Z, // sand slab values
+ int row, // target row
+ int col, // target col
+ int* row_max, // max slope neighbor row
+ int* col_max, // max slope neighbor col
+ double* dh_max) // max slope
+{
+
+ int n; // 1d neighbor idx
+ int i, j; // 2d neighbor idx
+ double dh; // slope
+
+ // relative neighbor coordinates and distances
+ int drow[] = {1, 1, 1, 0, -1, -1, -1, 0};
+ int dcol[] = {-1, 0, 1, 1, 1, 0, -1, -1};
+ double dx[] = {1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0};
+
+ // matrix dimensions
+ int nrows = Z->size1;
+ int ncols = Z->size2;
+
+ // check the 8 surrounding neighbor cells
+ for (n = 0; n<8; n++) {
+
+ // absolute neighbor cell idx
+ i = row + drow[n];
+ j = col + dcol[n];
+
+ // correct for periodic boundaries
+ if (i < 0)
+ i += nrows;
+ if (i >= nrows)
+ i -= nrows;
+ if (j < 0)
+ j += ncols;
+ if (j >= ncols)
+ j -= ncols;
+
+ // find slope
+ dh = (gsl_matrix_get(Z, i, j) - gsl_matrix_get(Z, row, col)) / dx[n];
+
+ // save position if it is the highest slope so far
+ if (dh > *dh_max) {
+ *dh_max = dh;
+ *row_max = i;
+ *col_max = j;
+ }
+ }
+}
+
+void find_max_slope_neighbor_depo(
+ gsl_matrix* Z, // sand slab values
+ int row, // target row
+ int col, // target col
+ int* row_max, // max slope neighbor row
+ int* col_max, // max slope neighbor col
+ double* dh_max) // max slope
+{
+
+ int n; // 1d neighbor idx
+ int i, j; // 2d neighbor idx
+ double dh; // slope
+
+ // relative neighbor coordinates and distances
+ int drow[] = {1, 1, 1, 0, -1, -1, -1, 0};
+ int dcol[] = {-1, 0, 1, 1, 1, 0, -1, -1};
+ double dx[] = {1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0};
+
+ // matrix dimensions
+ int nrows = Z->size1;
+ int ncols = Z->size2;
+
+ // check the 8 surrounding neighbor cells
+ for (n = 0; n<8; n++) {
+
+ // absolute neighbor cell idx
+ i = row + drow[n];
+ j = col + dcol[n];
+
+ // correct for periodic boundaries
+ if (i < 0)
+ i += nrows;
+ if (i >= nrows)
+ i -= nrows;
+ if (j < 0)
+ j += ncols;
+ if (j >= ncols)
+ j -= ncols;
+
+ // find slope
+ dh = (gsl_matrix_get(Z, row, col) - gsl_matrix_get(Z, i, j)) / dx[n];
+
+ // save position if it is the highest slope so far
+ if (dh > *dh_max) {
+ *dh_max = dh;
+ *row_max = i;
+ *col_max = j;
+ }
+ }
+
+}
+
+
+// Check and perform avalanche into cell if slope exceeds limit
+void avalanche_erosion(
+ gsl_matrix* Z, // sand slab values
+ int row, // target row
+ int col) // target col
+{
+ // find the neighbor with the max. height difference and the slope value
+ double dh_max = 0.0; // max slope
+ int row_max, col_max;
+ find_max_slope_neighbor_ero(Z, row, col, &row_max, &col_max, &dh_max);
+
+ // Perform avalanche along max slope neighbor
+ double threshold = 2.0; // avalanche threshold
+ if (dh_max > threshold) {
+ deposit(Z, row, col); // put sand in target cell
+ erode(Z, row_max, col_max); // remove sand from max slope neighbor cell
+
+ // Recursion
+ avalanche_erosion(Z, row, col);
+ }
+}
+
+// Check and perform avalanche away from cell if slope exceeds limit
+void avalanche_deposition(
+ gsl_matrix* Z, // sand slab values
+ int row, // target row
+ int col) // target col
+{
+ // find the neighbor with the max. height difference and the slope value
+ double dh_max = 0.0; // max slope
+ int row_max, col_max;
+ find_max_slope_neighbor_depo(Z, row, col, &row_max, &col_max, &dh_max);
+
+ // Perform avalanche along max slope neighbor
+ double threshold = 2.0; // avalanche threshold
+ if (dh_max > threshold) {
+ erode(Z, row, col);
+ deposit(Z, row_max, col_max);
+
+ // Recursion
+ avalanche_deposition(Z, row, col);
+ }
+}
+
+// Move a sand unit and deposit it where it fits
+void move_and_deposit(
+ gsl_matrix* Z, // matrix with sand slabs
+ double p_ns, // likelihood of deposition in no-sand cell
+ double p_s, // likelihood of deposition in sand cell
+ double d_l, // transport distance in no. of cells
+ int* row, // current cell row
+ int* col, // current cell col
+ int* deposited) // deposited? 1=yes, 0=no
+{
+ // move sand slab in wind direction
+ *col += d_l;
+ int ncols = Z->size1;
+ if (*col >= ncols)
+ *col -= ncols;
+
+ // is the target in the shade?
+ // 1=yes, 0=no
+ int shade = 0;
+ shade = inshade(Z, *row, *col);
+ if (shade > 0) {
+ deposit(Z, *row, *col);
+ avalanche_deposition(Z, *row, *col); // check for avalanches
+ *deposited = 1; // sand deposited, stop moving it
+ }
+
+ // if not in the shade, check against deposition probabilities
+ else {
+
+ // sand in cell
+ if (gsl_matrix_get(Z, *row, *col) > 0.) {
+ if (rnd() <= p_s) { // deposit in cell with sand
+ deposit(Z, *row, *col);
+ avalanche_deposition(Z, *row, *col); // check for avalanches
+ *deposited = 1; // sand deposited, stop moving it
+ }
+ } else { // no sand in cell
+ if (rnd() <= p_ns) { // deposit in cell without sand
+ deposit(Z, *row, *col);
+ avalanche_deposition(Z, *row, *col); // check for avalanches
+ *deposited = 1; // sand deposited, stop moving it
+ }
+ }
+ }
+}
+
+// Perform a single Werner iteration
+void werner_iter(
+ gsl_matrix* Z, // matrix with sand slabs
+ int d_l, // wind transport distance
+ double p_ns, // likelihood of deposition in sand-free cell
+ double p_s) // likelihood of deposition in sand-covered cell
+{
+ // Evaluate N=nx*ny cells in random order
+ int row, col, deposited;
+ int nrows = Z->size1; // row major
+ int ncols = Z->size2;
+ long int n;
+ long int N = nrows*ncols;
+ double cellval;
+
+
+ // Perform cycle
+ for (n=0; n<N; n++) {
+
+ // random cell idx
+ row = rand()%nrows;
+ col = rand()%ncols;
+
+ // If the cell has sand in it (val>0), and is not in the shadow zone
+ cellval = gsl_matrix_get(Z, row, col);
+ if ((cellval > 0.) && (inshade(Z, row, col) == 0)) {
+
+ // erode
+ erode(Z, row, col);
+
+ // check for avalanche
+ avalanche_erosion(Z, row, col);
+
+ // move eroded sand and deposit where it fits
+ deposited = 0;
+ while (deposited == 0)
+ move_and_deposit(Z, p_ns, p_s, d_l, &row, &col, &deposited);
+
+ }
+ }
+}
+
+// Loop system through time
+void werner_loop(
+ gsl_matrix* Z, // matrix with sand slabs
+ long int N_c, // number of iterations
+ int d_l, // wind transport distance
+ double p_ns, // likelihood of deposition in sand-free cell
+ double p_s) // likelihood of deposition in sand-covered cell
+{
+ long int t;
+ for (t = 0; t<N_c; t++)
+ werner_iter(Z, d_l, p_ns, p_s);
+}
+
+
(DIR) diff --git a/werner.h b/werner.h
t@@ -0,0 +1,23 @@
+#ifndef WERNER_H_
+#define WERNER_H_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <gsl/gsl_matrix.h>
+
+void init_rnd_matrix(gsl_matrix* M, double low, double high);
+void print_matrix(gsl_matrix* M);
+void werner_iter(
+ gsl_matrix* Z, // matrix with sand slabs
+ int d_l, // wind transport distance
+ double p_ns, // likelihood of deposition in sand-free cell
+ double p_s); // likelihood of deposition in sand-covered cell
+void werner_loop(
+ gsl_matrix* Z, // matrix with sand slabs
+ long int N_c, // number of iterations
+ int d_l, // wind transport distance
+ double p_ns, // likelihood of deposition in sand-free cell
+ double p_s); // likelihood of deposition in sand-covered cell
+
+
+#endif
(DIR) diff --git a/wernerparams.h b/wernerparams.h
t@@ -0,0 +1,11 @@
+#ifndef WERNERPARAMS_H_
+#define WERNERPARAMS_H_
+
+// set simulation parameters
+const int rows = 400;
+const int cols = 400;
+const int d_l = 5; // wind transport distance (left to right)
+const double p_ns = 0.4; // likelihood of deposition in sand free cells
+const double p_s = 0.6; // likelihood of deposition in sand covered cells
+
+#endif