twerner.c - werner - cellular automata simulation of wind-driven sand transport
HTML git clone git://src.adamsgaard.dk/werner
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
twerner.c (9111B)
---
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <gsl/gsl_matrix.h>
4
5 // Return a random number in [0;1[
6 double rnd(void)
7 {
8 return (double)rand()/RAND_MAX;
9 }
10
11 // Initialize matrix M with random ints from low to high
12 void init_rnd_matrix(gsl_matrix* M, double low, double high)
13 {
14 int row, col;
15 for (row = 0; row < M->size1; row++)
16 for (col = 0; col < M->size2; col++)
17 gsl_matrix_set(M, row, col,
18 (int)(rnd() * (high-low) + low));
19 }
20
21 // Write matrix to stdout
22 void print_matrix(gsl_matrix* M)
23 {
24 int row, col;
25 double val;
26 for (row = 0; row < M->size1; row++) {
27 for (col = 0; col < M->size2; col++) {
28 val = gsl_matrix_get(M, row, col);
29 printf("%f\t", val);
30 }
31 puts(""); // newline
32 }
33 }
34
35 // Add a slab of sand from the target cell
36 void deposit(gsl_matrix* Z, int row, int col)
37 {
38 gsl_matrix_set(Z, row, col, gsl_matrix_get(Z, row, col) + 1.0);
39 }
40
41 // Remove a slab of sand from the target cell
42 void erode(gsl_matrix* Z, int row, int col)
43 {
44 gsl_matrix_set(Z, row, col, gsl_matrix_get(Z, row, col) - 1.0);
45 }
46
47 // Find out if cell is in shade (1) or not (0)
48 int inshade(
49 gsl_matrix* Z, // matrix with sand slabs
50 int row, // row idx of interest
51 int col) // col idx of interest
52 {
53 int js;
54 int i = 1;
55 int shade = 0;
56 int nx = Z->size2;
57
58 while ((i <= nx/4) && (shade == 0)) {
59 js = col - i;
60 if (js < 0) // periodic boundary
61 js += nx;
62 if (gsl_matrix_get(Z, row, js) >= gsl_matrix_get(Z, row, col) + i)
63 shade = 1;
64 i++;
65 }
66 return shade;
67 }
68
69 void find_max_slope_neighbor_ero(
70 gsl_matrix* Z, // sand slab values
71 int row, // target row
72 int col, // target col
73 int* row_max, // max slope neighbor row
74 int* col_max, // max slope neighbor col
75 double* dh_max) // max slope
76 {
77
78 int n; // 1d neighbor idx
79 int i, j; // 2d neighbor idx
80 double dh; // slope
81
82 // relative neighbor coordinates and distances
83 int drow[] = {1, 1, 1, 0, -1, -1, -1, 0};
84 int dcol[] = {-1, 0, 1, 1, 1, 0, -1, -1};
85 double dx[] = {1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0};
86
87 // matrix dimensions
88 int nrows = Z->size1;
89 int ncols = Z->size2;
90
91 // check the 8 surrounding neighbor cells
92 for (n = 0; n<8; n++) {
93
94 // absolute neighbor cell idx
95 i = row + drow[n];
96 j = col + dcol[n];
97
98 // correct for periodic boundaries
99 if (i < 0)
100 i += nrows;
101 if (i >= nrows)
102 i -= nrows;
103 if (j < 0)
104 j += ncols;
105 if (j >= ncols)
106 j -= ncols;
107
108 // find slope
109 dh = (gsl_matrix_get(Z, i, j) - gsl_matrix_get(Z, row, col)) / dx[n];
110
111 // save position if it is the highest slope so far
112 if (dh > *dh_max) {
113 *dh_max = dh;
114 *row_max = i;
115 *col_max = j;
116 }
117 }
118 }
119
120 void find_max_slope_neighbor_depo(
121 gsl_matrix* Z, // sand slab values
122 int row, // target row
123 int col, // target col
124 int* row_max, // max slope neighbor row
125 int* col_max, // max slope neighbor col
126 double* dh_max) // max slope
127 {
128
129 int n; // 1d neighbor idx
130 int i, j; // 2d neighbor idx
131 double dh; // slope
132
133 // relative neighbor coordinates and distances
134 int drow[] = {1, 1, 1, 0, -1, -1, -1, 0};
135 int dcol[] = {-1, 0, 1, 1, 1, 0, -1, -1};
136 double dx[] = {1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0, 1.4142, 1.0};
137
138 // matrix dimensions
139 int nrows = Z->size1;
140 int ncols = Z->size2;
141
142 // check the 8 surrounding neighbor cells
143 for (n = 0; n<8; n++) {
144
145 // absolute neighbor cell idx
146 i = row + drow[n];
147 j = col + dcol[n];
148
149 // correct for periodic boundaries
150 if (i < 0)
151 i += nrows;
152 if (i >= nrows)
153 i -= nrows;
154 if (j < 0)
155 j += ncols;
156 if (j >= ncols)
157 j -= ncols;
158
159 // find slope
160 dh = (gsl_matrix_get(Z, row, col) - gsl_matrix_get(Z, i, j)) / dx[n];
161
162 // save position if it is the highest slope so far
163 if (dh > *dh_max) {
164 *dh_max = dh;
165 *row_max = i;
166 *col_max = j;
167 }
168 }
169 }
170
171 // Check and perform avalanche into cell if slope exceeds limit
172 void avalanche_erosion(
173 gsl_matrix* Z, // sand slab values
174 int row, // target row
175 int col) // target col
176 {
177 // find the neighbor with the max. height difference and the slope value
178 double dh_max = 0.0; // max slope
179 int row_max, col_max;
180 find_max_slope_neighbor_ero(Z, row, col, &row_max, &col_max, &dh_max);
181
182 // Perform avalanche along max slope neighbor
183 double threshold = 2.0; // avalanche threshold
184 if (dh_max > threshold) {
185 deposit(Z, row, col); // put sand in target cell
186 erode(Z, row_max, col_max); // remove sand from max slope neighbor cell
187
188 // Recursion
189 avalanche_erosion(Z, row, col);
190 }
191 }
192
193 // Check and perform avalanche away from cell if slope exceeds limit
194 void avalanche_deposition(
195 gsl_matrix* Z, // sand slab values
196 int row, // target row
197 int col) // target col
198 {
199 // find the neighbor with the max. height difference and the slope value
200 double dh_max = 0.0; // max slope
201 int row_max, col_max;
202 find_max_slope_neighbor_depo(Z, row, col, &row_max, &col_max, &dh_max);
203
204 // Perform avalanche along max slope neighbor
205 double threshold = 2.0; // avalanche threshold
206 if (dh_max > threshold) {
207 erode(Z, row, col);
208 deposit(Z, row_max, col_max);
209
210 // Recursion
211 avalanche_deposition(Z, row, col);
212 }
213 }
214
215 // Move a sand unit and deposit it where it fits
216 void move_and_deposit(
217 gsl_matrix* Z, // matrix with sand slabs
218 double p_ns, // likelihood of deposition in no-sand cell
219 double p_s, // likelihood of deposition in sand cell
220 double d_l, // transport distance in no. of cells
221 int* row, // current cell row
222 int* col, // current cell col
223 int* deposited) // deposited? 1=yes, 0=no
224 {
225 // move sand slab in wind direction
226 *col += d_l;
227 int ncols = Z->size1;
228 if (*col >= ncols)
229 *col -= ncols;
230
231 // is the target in the shade?
232 // 1=yes, 0=no
233 int shade = 0;
234 shade = inshade(Z, *row, *col);
235 if (shade > 0) {
236 deposit(Z, *row, *col);
237 avalanche_deposition(Z, *row, *col); // check for avalanches
238 *deposited = 1; // sand deposited, stop moving it
239 }
240
241 // if not in the shade, check against deposition probabilities
242 else {
243
244 // sand in cell
245 if (gsl_matrix_get(Z, *row, *col) > 0.) {
246 if (rnd() <= p_s) { // deposit in cell with sand
247 deposit(Z, *row, *col);
248 avalanche_deposition(Z, *row, *col); // check for avalanches
249 *deposited = 1; // sand deposited, stop moving it
250 }
251 } else { // no sand in cell
252 if (rnd() <= p_ns) { // deposit in cell without sand
253 deposit(Z, *row, *col);
254 avalanche_deposition(Z, *row, *col); // check for avalanches
255 *deposited = 1; // sand deposited, stop moving it
256 }
257 }
258 }
259 }
260
261 // Perform a single Werner iteration
262 void werner_iter(
263 gsl_matrix* Z, // matrix with sand slabs
264 int d_l, // wind transport distance
265 double p_ns, // likelihood of deposition in sand-free cell
266 double p_s) // likelihood of deposition in sand-covered cell
267 {
268 // Evaluate N=nx*ny cells in random order
269 int row, col, deposited;
270 int nrows = Z->size1; // row major
271 int ncols = Z->size2;
272 long int n;
273 long int N = nrows*ncols;
274 double cellval;
275
276
277 // Perform cycle
278 for (n=0; n<N; n++) {
279
280 // random cell idx
281 row = rand()%nrows;
282 col = rand()%ncols;
283
284 // If the cell has sand in it (val>0), and is not in the shadow zone
285 cellval = gsl_matrix_get(Z, row, col);
286 if ((cellval > 0.) && (inshade(Z, row, col) == 0)) {
287
288 // erode
289 erode(Z, row, col);
290
291 // check for avalanche
292 avalanche_erosion(Z, row, col);
293
294 // move eroded sand and deposit where it fits
295 deposited = 0;
296 while (deposited == 0)
297 move_and_deposit(Z, p_ns, p_s, d_l, &row, &col, &deposited);
298
299 }
300 }
301 }
302
303 // Loop system through time
304 void werner_loop(
305 gsl_matrix* Z, // matrix with sand slabs
306 long int N_c, // number of iterations
307 int d_l, // wind transport distance
308 double p_ns, // likelihood of deposition in sand-free cell
309 double p_s) // likelihood of deposition in sand-covered cell
310 {
311 long int t;
312 for (t = 0; t<N_c; t++)
313 werner_iter(Z, d_l, p_ns, p_s);
314 }