Context Navigation

source: anuga_core/source/anuga/shallow_water/shallow_water_ext.c @ 4815

Last change on this file since 4815 was 4815, checked in by ole, 16 years ago
Work towards Froude_number adjusted, balanced limiters. This is currently only in place for tight_slope_limiters == 1.
File size: 70.7 KB

Line
1	// Python - C extension module for shallow_water.py
2	//
3	// To compile (Python2.3):
4	// gcc -c domain_ext.c -I/usr/include/python2.3 -o domain_ext.o -Wall -O
5	// gcc -shared domain_ext.o -o domain_ext.so
6	//
7	// or use python compile.py
8	//
9	// See the module shallow_water.py
10	//
11	//
12	// Ole Nielsen, GA 2004
13
14
15	#include "Python.h"
16	#include "Numeric/arrayobject.h"
17	#include "math.h"
18	#include <stdio.h>
19
20	// Shared code snippets
21	#include "util_ext.h"
22
23
24	const double pi = 3.14159265358979;
25
26	// Computational function for rotation
27	int _rotate(double *q, double n1, double n2) {
28	/*Rotate the momentum component q (q[1], q[2])
29	from x,y coordinates to coordinates based on normal vector (n1, n2).
30
31	Result is returned in array 3x1 r
32	To rotate in opposite direction, call rotate with (q, n1, -n2)
33
34	Contents of q are changed by this function */
35
36
37	double q1, q2;
38
39	// Shorthands
40	q1 = q[1]; // uh momentum
41	q2 = q[2]; // vh momentum
42
43	// Rotate
44	q[1] = n1q1 + n2q2;
45	q[2] = -n2q1 + n1q2;
46
47	return 0;
48	}
49
50	int find_qmin_and_qmax(double dq0, double dq1, double dq2,
51	double qmin, double qmax){
52	// Considering the centroid of an FV triangle and the vertices of its
53	// auxiliary triangle, find
54	// qmin=min(q)-qc and qmax=max(q)-qc,
55	// where min(q) and max(q) are respectively min and max over the
56	// four values (at the centroid of the FV triangle and the auxiliary
57	// triangle vertices),
58	// and qc is the centroid
59	// dq0=q(vertex0)-q(centroid of FV triangle)
60	// dq1=q(vertex1)-q(vertex0)
61	// dq2=q(vertex2)-q(vertex0)
62
63	if (dq0>=0.0){
64	if (dq1>=dq2){
65	if (dq1>=0.0)
66	*qmax=dq0+dq1;
67	else
68	*qmax=dq0;
69	if ((*qmin=dq0+dq2)<0)
70	;// qmin is already set to correct value
71	else
72	*qmin=0.0;
73	}
74	else{// dq1<dq2
75	if (dq2>0)
76	*qmax=dq0+dq2;
77	else
78	*qmax=dq0;
79	if ((*qmin=dq0+dq1)<0)
80	;// qmin is the correct value
81	else
82	*qmin=0.0;
83	}
84	}
85	else{//dq0<0
86	if (dq1<=dq2){
87	if (dq1<0.0)
88	*qmin=dq0+dq1;
89	else
90	*qmin=dq0;
91	if ((*qmax=dq0+dq2)>0.0)
92	;// qmax is already set to the correct value
93	else
94	*qmax=0.0;
95	}
96	else{// dq1>dq2
97	if (dq2<0.0)
98	*qmin=dq0+dq2;
99	else
100	*qmin=dq0;
101	if ((*qmax=dq0+dq1)>0.0)
102	;// qmax is already set to the correct value
103	else
104	*qmax=0.0;
105	}
106	}
107	return 0;
108	}
109
110	int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
111	// Given provisional jumps dqv from the FV triangle centroid to its
112	// vertices and jumps qmin (qmax) between the centroid of the FV
113	// triangle and the minimum (maximum) of the values at the centroid of
114	// the FV triangle and the auxiliary triangle vertices,
115	// calculate a multiplicative factor phi by which the provisional
116	// vertex jumps are to be limited
117
118	int i;
119	double r=1000.0, r0=1.0, phi=1.0;
120	static double TINY = 1.0e-100; // to avoid machine accuracy problems.
121	// FIXME: Perhaps use the epsilon used elsewhere.
122
123	// Any provisional jump with magnitude < TINY does not contribute to
124	// the limiting process.
125	for (i=0;i<3;i++){
126	if (dqv[i]<-TINY)
127	r0=qmin/dqv[i];
128
129	if (dqv[i]>TINY)
130	r0=qmax/dqv[i];
131
132	r=min(r0,r);
133	}
134
135	phi=min(r*beta_w,1.0);
136	for (i=0;i<3;i++)
137	dqv[i]=dqv[i]*phi;
138
139	return 0;
140	}
141
142
143	double compute_froude_number(double uh,
144	double h,
145	double g,
146	double epsilon) {
147
148	// Compute Froude number; v/sqrt(gh)
149
150	double froude_number;
151
152	//Compute Froude number (stability diagnostics)
153	if (h > epsilon) {
154	froude_number = uh/sqrt(g*h)/h;
155	} else {
156	froude_number = 0.0;
157	// FIXME (Ole): What should it be when dry??
158	}
159
160	return froude_number;
161	}
162
163
164
165	// Function to obtain speed from momentum and depth.
166	// This is used by flux functions
167	// Input parameters uh and h may be modified by this function.
168	double _compute_speed(double *uh,
169	double *h,
170	double epsilon,
171	double h0) {
172
173	double u;
174
175
176	if (*h < epsilon) {
177	*h = 0.0; //Could have been negative
178	u = 0.0;
179	} else {
180	u = uh/(h + h0/ *h);
181	}
182
183
184	// Adjust momentum to be consistent with speed
185	uh = u *h;
186
187	return u;
188	}
189
190	// Innermost flux function (using stage w=z+h)
191	int _flux_function_central(double q_left, double q_right,
192	double z_left, double z_right,
193	double n1, double n2,
194	double epsilon, double H0, double g,
195	double edgeflux, double max_speed) {
196
197	/*Compute fluxes between volumes for the shallow water wave equation
198	cast in terms of the 'stage', w = h+z using
199	the 'central scheme' as described in
200
201	Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
202	Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
203	Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
204
205	The implemented formula is given in equation (3.15) on page 714
206	*/
207
208	int i;
209
210	double w_left, h_left, uh_left, vh_left, u_left;
211	double w_right, h_right, uh_right, vh_right, u_right;
212	double v_left, v_right;
213	double s_min, s_max, soundspeed_left, soundspeed_right;
214	double denom, z;
215
216	// Workspace (allocate once, use many)
217	static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
218
219	double h0 = H0*H0; // This ensures a good balance when h approaches H0.
220	// But evidence suggests that h0 can be as little as
221	// epsilon!
222
223	// Copy conserved quantities to protect from modification
224	for (i=0; i<3; i++) {
225	q_left_rotated[i] = q_left[i];
226	q_right_rotated[i] = q_right[i];
227	}
228
229	// Align x- and y-momentum with x-axis
230	_rotate(q_left_rotated, n1, n2);
231	_rotate(q_right_rotated, n1, n2);
232
233	z = (z_left+z_right)/2; // Average elevation values
234
235	// Compute speeds in x-direction
236	w_left = q_left_rotated[0];
237	h_left = w_left-z;
238	uh_left = q_left_rotated[1];
239	u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
240
241	w_right = q_right_rotated[0];
242	h_right = w_right-z;
243	uh_right = q_right_rotated[1];
244	u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
245
246	// Momentum in y-direction
247	vh_left = q_left_rotated[2];
248	vh_right = q_right_rotated[2];
249
250	// Limit y-momentum if necessary
251	v_left = _compute_speed(&vh_left, &h_left, epsilon, h0);
252	v_right = _compute_speed(&vh_right, &h_right, epsilon, h0);
253
254	// Maximal and minimal wave speeds
255	soundspeed_left = sqrt(g*h_left);
256	soundspeed_right = sqrt(g*h_right);
257
258	s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
259	if (s_max < 0.0) s_max = 0.0;
260
261	s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
262	if (s_min > 0.0) s_min = 0.0;
263
264
265	// Flux formulas
266	flux_left[0] = u_left*h_left;
267	flux_left[1] = u_leftuh_left + 0.5gh_lefth_left;
268	flux_left[2] = u_left*vh_left;
269
270	flux_right[0] = u_right*h_right;
271	flux_right[1] = u_rightuh_right + 0.5gh_righth_right;
272	flux_right[2] = u_right*vh_right;
273
274
275	// Flux computation
276	denom = s_max-s_min;
277	if (denom < epsilon) { // FIXME (Ole): Try using H0 here
278	for (i=0; i<3; i++) edgeflux[i] = 0.0;
279	*max_speed = 0.0;
280	} else {
281	for (i=0; i<3; i++) {
282	edgeflux[i] = s_maxflux_left[i] - s_minflux_right[i];
283	edgeflux[i] += s_maxs_min(q_right_rotated[i]-q_left_rotated[i]);
284	edgeflux[i] /= denom;
285	}
286
287	// Maximal wavespeed
288	*max_speed = max(fabs(s_max), fabs(s_min));
289
290	// Rotate back
291	_rotate(edgeflux, n1, -n2);
292	}
293
294	return 0;
295	}
296
297	double erfcc(double x){
298	double z,t,result;
299
300	z=fabs(x);
301	t=1.0/(1.0+0.5*z);
302	result=texp(-zz-1.26551223+t(1.00002368+t(.37409196+
303	t(.09678418+t(-.18628806+t(.27886807+t(-1.13520398+
304	t(1.48851587+t(-.82215223+t*.17087277)))))))));
305	if (x < 0.0) result = 2.0-result;
306
307	return result;
308	}
309
310
311
312	// Computational function for flux computation (using stage w=z+h)
313	int flux_function_kinetic(double q_left, double q_right,
314	double z_left, double z_right,
315	double n1, double n2,
316	double epsilon, double H0, double g,
317	double edgeflux, double max_speed) {
318
319	/*Compute fluxes between volumes for the shallow water wave equation
320	cast in terms of the 'stage', w = h+z using
321	the 'central scheme' as described in
322
323	Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647.
324	*/
325
326	int i;
327
328	double w_left, h_left, uh_left, vh_left, u_left, F_left;
329	double w_right, h_right, uh_right, vh_right, u_right, F_right;
330	double s_min, s_max, soundspeed_left, soundspeed_right;
331	double z;
332	double q_left_rotated[3], q_right_rotated[3];
333
334	double h0 = H0*H0; //This ensures a good balance when h approaches H0.
335
336	//Copy conserved quantities to protect from modification
337	for (i=0; i<3; i++) {
338	q_left_rotated[i] = q_left[i];
339	q_right_rotated[i] = q_right[i];
340	}
341
342	//Align x- and y-momentum with x-axis
343	_rotate(q_left_rotated, n1, n2);
344	_rotate(q_right_rotated, n1, n2);
345
346	z = (z_left+z_right)/2; //Take average of field values
347
348	//Compute speeds in x-direction
349	w_left = q_left_rotated[0];
350	h_left = w_left-z;
351	uh_left = q_left_rotated[1];
352	u_left =_compute_speed(&uh_left, &h_left, epsilon, h0);
353
354	w_right = q_right_rotated[0];
355	h_right = w_right-z;
356	uh_right = q_right_rotated[1];
357	u_right =_compute_speed(&uh_right, &h_right, epsilon, h0);
358
359
360	//Momentum in y-direction
361	vh_left = q_left_rotated[2];
362	vh_right = q_right_rotated[2];
363
364
365	//Maximal and minimal wave speeds
366	soundspeed_left = sqrt(g*h_left);
367	soundspeed_right = sqrt(g*h_right);
368
369	s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
370	if (s_max < 0.0) s_max = 0.0;
371
372	s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
373	if (s_min > 0.0) s_min = 0.0;
374
375
376	F_left = 0.0;
377	F_right = 0.0;
378	if (h_left > 0.0) F_left = u_left/sqrt(g*h_left);
379	if (h_right > 0.0) F_right = u_right/sqrt(g*h_right);
380
381	for (i=0; i<3; i++) edgeflux[i] = 0.0;
382	*max_speed = 0.0;
383
384	edgeflux[0] = h_leftu_left/2.0erfcc(-F_left) + \
385	h_leftsqrt(gh_left)/2.0/sqrt(pi)exp(-(F_leftF_left)) + \
386	h_rightu_right/2.0erfcc(F_right) - \
387	h_rightsqrt(gh_right)/2.0/sqrt(pi)exp(-(F_rightF_right));
388
389	edgeflux[1] = (h_leftu_leftu_left + g/2.0h_lefth_left)/2.0*erfcc(-F_left) + \
390	u_lefth_leftsqrt(gh_left)/2.0/sqrt(pi)exp(-(F_left*F_left)) + \
391	(h_rightu_rightu_right + g/2.0h_righth_right)/2.0*erfcc(F_right) - \
392	u_righth_rightsqrt(gh_right)/2.0/sqrt(pi)exp(-(F_right*F_right));
393
394	edgeflux[2] = vh_leftu_left/2.0erfcc(-F_left) + \
395	vh_leftsqrt(gh_left)/2.0/sqrt(pi)exp(-(F_leftF_left)) + \
396	vh_rightu_right/2.0erfcc(F_right) - \
397	vh_rightsqrt(gh_right)/2.0/sqrt(pi)exp(-(F_rightF_right));
398
399	//Maximal wavespeed
400	*max_speed = max(fabs(s_max), fabs(s_min));
401
402	//Rotate back
403	_rotate(edgeflux, n1, -n2);
404
405	return 0;
406	}
407
408
409
410
411	void _manning_friction(double g, double eps, int N,
412	double* w, double* z,
413	double* uh, double* vh,
414	double* eta, double* xmom, double* ymom) {
415
416	int k;
417	double S, h;
418
419	for (k=0; k<N; k++) {
420	if (eta[k] > eps) {
421	h = w[k]-z[k];
422	if (h >= eps) {
423	S = -g * eta[k]eta[k] sqrt((uh[k]uh[k] + vh[k]vh[k]));
424	S /= pow(h, 7.0/3); //Expensive (on Ole's home computer)
425	//S /= exp(7.0/3.0*log(h)); //seems to save about 15% over manning_friction
426	//S /= hh(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
427
428
429	//Update momentum
430	xmom[k] += S*uh[k];
431	ymom[k] += S*vh[k];
432	}
433	}
434	}
435	}
436
437
438	/*
439	void _manning_friction_explicit(double g, double eps, int N,
440	double* w, double* z,
441	double* uh, double* vh,
442	double* eta, double* xmom, double* ymom) {
443
444	int k;
445	double S, h;
446
447	for (k=0; k<N; k++) {
448	if (eta[k] > eps) {
449	h = w[k]-z[k];
450	if (h >= eps) {
451	S = -g * eta[k]eta[k] sqrt((uh[k]uh[k] + vh[k]vh[k]));
452	S /= pow(h, 7.0/3); //Expensive (on Ole's home computer)
453	//S /= exp(7.0/3.0*log(h)); //seems to save about 15% over manning_friction
454	//S /= hh(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
455
456
457	//Update momentum
458	xmom[k] += S*uh[k];
459	ymom[k] += S*vh[k];
460	}
461	}
462	}
463	}
464	*/
465
466	int _balance_deep_and_shallow(int N,
467	double beta_h,
468	double* wc,
469	double* zc,
470	double* wv,
471	double* zv,
472	double* hvbar, // Retire this
473	double* xmomc,
474	double* ymomc,
475	double* xmomv,
476	double* ymomv,
477	double H0,
478	int tight_slope_limiters,
479	double alpha_balance) {
480
481	int k, k3, i, excessive_froude_number=0;
482	double dz, hmin, alpha, h_diff, hc_k;
483	double epsilon = 1.0e-6; // FIXME: Temporary measure
484	double g = 9.81; // FIXME: Temporary measure
485	double hv[3], h; // Depths at vertices
486	double Fx, Fy; // Froude numbers
487
488	// Compute linear combination between w-limited stages and
489	// h-limited stages close to the bed elevation.
490
491	for (k=0; k<N; k++) {
492	// Compute maximal variation in bed elevation
493	// This quantitiy is
494	// dz = max_i abs(z_i - z_c)
495	// and it is independent of dimension
496	// In the 1d case zc = (z0+z1)/2
497	// In the 2d case zc = (z0+z1+z2)/3
498
499	k3 = 3*k;
500	hc_k = wc[k] - zc[k]; // Centroid value at triangle k
501
502	dz = 0.0;
503	if (tight_slope_limiters == 0) {
504	// FIXME: Try with this one precomputed
505	for (i=0; i<3; i++) {
506	dz = max(dz, fabs(zv[k3+i]-zc[k]));
507	}
508	}
509
510	// Calculate depth at vertices (possibly negative here!)
511	hv[0] = wv[k3] - zv[k3];
512	hv[1] = wv[k3+1] - zv[k3+1];
513	hv[2] = wv[k3+2] - zv[k3+2];
514
515	// Calculate minimal depth across all three vertices
516	hmin = min(hv[0], min(hv[1], hv[2]));
517
518
519	// Create alpha in [0,1], where alpha==0 means using the h-limited
520	// stage and alpha==1 means using the w-limited stage as
521	// computed by the gradient limiter (both 1st or 2nd order)
522	if (tight_slope_limiters == 0) {
523	// If hmin > dz/alpha_balance then alpha = 1 and the bed will have no
524	// effect
525	// If hmin < 0 then alpha = 0 reverting to constant height above bed.
526	// The parameter alpha_balance==2 by default
527
528
529	if (dz > 0.0) {
530	alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );
531	} else {
532	alpha = 1.0; // Flat bed
533	}
534	//printf("Using old style limiter\n");
535
536	} else {
537
538	// Tight Slope Limiter (2007)
539
540	// Make alpha as large as possible but still ensure that
541	// final depth is positive
542
543	if (hmin < H0) {
544	alpha = 1.0;
545	for (i=0; i<3; i++) {
546
547	// FIXME (Ole): Simplify when (if) hvbar gets retired
548	if (beta_h > epsilon) {
549	h_diff = hvbar[k3+i] - hv[i];
550	} else {
551	h_diff = hc_k - hv[i];
552	}
553
554	if (h_diff <= 0) {
555	// Deep water triangle is further away from bed than
556	// shallow water (hbar < h). Any alpha will do
557
558	} else {
559	// Denominator is positive which means that we need some of the
560	// h-limited stage.
561
562	// FIXME (Ole): Simplify when (if) hvbar gets retired
563	if (beta_h > epsilon) {
564	alpha = min(alpha, (hvbar[k3+i] - H0)/h_diff);
565	} else {
566	alpha = min(alpha, (hc_k - H0)/h_diff);
567	}
568	}
569	}
570
571	// Ensure alpha in [0,1]
572	if (alpha>1.0) alpha=1.0;
573	if (alpha<0.0) alpha=0.0;
574
575	} else {
576	// Use w-limited stage exclusively
577	alpha = 1.0;
578	}
579	}
580
581
582
583	//printf("k=%d, hmin=%.2f, dz=%.2f, alpha=%.2f, alpha_balance=%.2f\n",
584	// k, hmin, dz, alpha, alpha_balance);
585
586	//printf("dz = %.3f, alpha = %.8f\n", dz, alpha);
587
588	// Let
589	//
590	// wvi be the w-limited stage (wvi = zvi + hvi)
591	// wvi- be the h-limited state (wvi- = zvi + hvi-)
592	//
593	//
594	// where i=0,1,2 denotes the vertex ids
595	//
596	// Weighted balance between w-limited and h-limited stage is
597	//
598	// wvi := (1-alpha)(zvi+hvi-) + alpha(zvi+hvi)
599	//
600	// It follows that the updated wvi is
601	// wvi := zvi + (1-alpha)hvi- + alphahvi
602	//
603	// Momentum is balanced between constant and limited
604
605
606	if (alpha < 1) {
607	for (i=0; i<3; i++) {
608	// FIXME (Ole): Simplify when (if) hvbar gets retired
609	if (beta_h > epsilon) {
610	wv[k3+i] = zv[k3+i] + (1-alpha)hvbar[k3+i] + alphahv[i];
611	} else {
612	wv[k3+i] = zv[k3+i] + (1-alpha)hc_k + alphahv[i];
613	}
614
615	// Update momentum as a linear combination of
616	// xmomc and ymomc (shallow) and momentum
617	// from extrapolator xmomv and ymomv (deep).
618	// FIXME (Ole): Is this really needed?
619	xmomv[k3+i] = (1-alpha)xmomc[k] + alphaxmomv[k3+i];
620	ymomv[k3+i] = (1-alpha)ymomc[k] + alphaymomv[k3+i];
621	}
622	}
623
624
625
626	if (tight_slope_limiters == 1) {
627
628	// Ensure that the Froude number is kept realistic at vertices
629	// FIXME (Ole): I think it could be used to adjust alpha down
630	// whenever Fr is too large. Possible make sure it doesn't deviate
631	// too much from the value at the centroid. I like this idea!
632
633	// FIXME (Ole): currently only used with tights_SL
634
635	excessive_froude_number=0;
636	for (i=0; i<3; i++) {
637
638	// Recalculate depth at vertex i
639	h = wv[k3+i] - zv[k3+i];
640
641	Fx = compute_froude_number(xmomv[k3+i], h, g, epsilon);
642	Fy = compute_froude_number(ymomv[k3+i], h, g, epsilon);
643
644	if ( (fabs(Fx) > 100.0) \|\| (fabs(Fy) > 100.0)) {
645	// FIXME: use max_froude - or base it on centroid value of F
646	// printf("Excessive Froude number detected: %f or %f\n", Fx, Fy);
647	excessive_froude_number=1;
648	}
649	}
650
651	if (excessive_froude_number) {
652
653	// printf("Adjusting momentum to first order.\n");
654	// Go back to first order (alpha = 0) for this triangle
655	for (i=0; i<3; i++) {
656	xmomv[k3+i] = xmomc[k];
657	ymomv[k3+i] = ymomc[k];
658
659	if (beta_h > epsilon) {
660	wv[k3+i] = zv[k3+i] + hvbar[k3+i];
661	} else {
662	wv[k3+i] = zv[k3+i] + hc_k;
663	}
664	}
665	}
666	}
667	}
668	return 0;
669	}
670
671
672
673	int _protect(int N,
674	double minimum_allowed_height,
675	double maximum_allowed_speed,
676	double epsilon,
677	double* wc,
678	double* zc,
679	double* xmomc,
680	double* ymomc) {
681
682	int k;
683	double hc;
684	double u, v, reduced_speed;
685
686
687	// Protect against initesimal and negative heights
688	if (maximum_allowed_speed < epsilon) {
689	for (k=0; k<N; k++) {
690	hc = wc[k] - zc[k];
691
692	if (hc < minimum_allowed_height) {
693
694	// Set momentum to zero and ensure h is non negative
695	xmomc[k] = 0.0;
696	ymomc[k] = 0.0;
697	if (hc <= 0.0) wc[k] = zc[k];
698	}
699	}
700	} else {
701
702	// Protect against initesimal and negative heights
703	for (k=0; k<N; k++) {
704	hc = wc[k] - zc[k];
705
706	if (hc < minimum_allowed_height) {
707
708	//New code: Adjust momentum to guarantee speeds are physical
709	// ensure h is non negative
710
711	if (hc <= 0.0) {
712	wc[k] = zc[k];
713	xmomc[k] = 0.0;
714	ymomc[k] = 0.0;
715	} else {
716	//Reduce excessive speeds derived from division by small hc
717	//FIXME (Ole): This may be unnecessary with new slope limiters
718	//in effect.
719
720	u = xmomc[k]/hc;
721	if (fabs(u) > maximum_allowed_speed) {
722	reduced_speed = maximum_allowed_speed * u/fabs(u);
723	//printf("Speed (u) has been reduced from %.3f to %.3f\n",
724	// u, reduced_speed);
725	xmomc[k] = reduced_speed * hc;
726	}
727
728	v = ymomc[k]/hc;
729	if (fabs(v) > maximum_allowed_speed) {
730	reduced_speed = maximum_allowed_speed * v/fabs(v);
731	//printf("Speed (v) has been reduced from %.3f to %.3f\n",
732	// v, reduced_speed);
733	ymomc[k] = reduced_speed * hc;
734	}
735	}
736	}
737	}
738	}
739	return 0;
740	}
741
742
743
744	int _assign_wind_field_values(int N,
745	double* xmom_update,
746	double* ymom_update,
747	double* s_vec,
748	double* phi_vec,
749	double cw) {
750
751	// Assign windfield values to momentum updates
752
753	int k;
754	double S, s, phi, u, v;
755
756	for (k=0; k<N; k++) {
757
758	s = s_vec[k];
759	phi = phi_vec[k];
760
761	//Convert to radians
762	phi = phi*pi/180;
763
764	//Compute velocity vector (u, v)
765	u = s*cos(phi);
766	v = s*sin(phi);
767
768	//Compute wind stress
769	S = cw * sqrt(uu + vv);
770	xmom_update[k] += S*u;
771	ymom_update[k] += S*v;
772	}
773	return 0;
774	}
775
776
777
778	///////////////////////////////////////////////////////////////////
779	// Gateways to Python
780
781
782
783	PyObject flux_function_central(PyObject self, PyObject *args) {
784	//
785	// Gateway to innermost flux function.
786	// This is only used by the unit tests as the c implementation is
787	// normally called by compute_fluxes in this module.
788
789
790	PyArrayObject normal, ql, qr, edgeflux;
791	double g, epsilon, max_speed, H0, zl, zr;
792
793	if (!PyArg_ParseTuple(args, "OOOddOddd",
794	&normal, &ql, &qr, &zl, &zr, &edgeflux,
795	&epsilon, &g, &H0)) {
796	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
797	return NULL;
798	}
799
800
801	_flux_function_central((double*) ql -> data,
802	(double*) qr -> data,
803	zl,
804	zr,
805	((double*) normal -> data)[0],
806	((double*) normal -> data)[1],
807	epsilon, H0, g,
808	(double*) edgeflux -> data,
809	&max_speed);
810
811	return Py_BuildValue("d", max_speed);
812	}
813
814
815
816	PyObject gravity(PyObject self, PyObject *args) {
817	//
818	// gravity(g, h, v, x, xmom, ymom)
819	//
820
821
822	PyArrayObject h, v, x, xmom, *ymom;
823	int k, N, k3, k6;
824	double g, avg_h, zx, zy;
825	double x0, y0, x1, y1, x2, y2, z0, z1, z2;
826	//double epsilon;
827
828	if (!PyArg_ParseTuple(args, "dOOOOO",
829	&g, &h, &v, &x,
830	&xmom, &ymom)) {
831	//&epsilon)) {
832	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
833	return NULL;
834	}
835
836	N = h -> dimensions[0];
837	for (k=0; k<N; k++) {
838	k3 = 3*k; // base index
839
840	// Get bathymetry
841	z0 = ((double*) v -> data)[k3 + 0];
842	z1 = ((double*) v -> data)[k3 + 1];
843	z2 = ((double*) v -> data)[k3 + 2];
844
845	// Optimise for flat bed
846	// Note (Ole): This didn't produce measurable speed up.
847	// Revisit later
848	//if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
849	// continue;
850	//}
851
852	// Get average depth from centroid values
853	avg_h = ((double *) h -> data)[k];
854
855	// Compute bed slope
856	k6 = 6*k; // base index
857
858	x0 = ((double*) x -> data)[k6 + 0];
859	y0 = ((double*) x -> data)[k6 + 1];
860	x1 = ((double*) x -> data)[k6 + 2];
861	y1 = ((double*) x -> data)[k6 + 3];
862	x2 = ((double*) x -> data)[k6 + 4];
863	y2 = ((double*) x -> data)[k6 + 5];
864
865
866	_gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
867
868	// Update momentum
869	((double) xmom -> data)[k] += -gzx*avg_h;
870	((double) ymom -> data)[k] += -gzy*avg_h;
871	}
872
873	return Py_BuildValue("");
874	}
875
876
877	PyObject manning_friction(PyObject self, PyObject *args) {
878	//
879	// manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
880	//
881
882
883	PyArrayObject w, z, uh, vh, eta, xmom, *ymom;
884	int N;
885	double g, eps;
886
887	if (!PyArg_ParseTuple(args, "ddOOOOOOO",
888	&g, &eps, &w, &z, &uh, &vh, &eta,
889	&xmom, &ymom)) {
890	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
891	return NULL;
892	}
893
894
895	N = w -> dimensions[0];
896	_manning_friction(g, eps, N,
897	(double*) w -> data,
898	(double*) z -> data,
899	(double*) uh -> data,
900	(double*) vh -> data,
901	(double*) eta -> data,
902	(double*) xmom -> data,
903	(double*) ymom -> data);
904
905	return Py_BuildValue("");
906	}
907
908
909	/*
910	PyObject manning_friction_explicit(PyObject self, PyObject *args) {
911	//
912	// manning_friction_explicit(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
913	//
914
915
916	PyArrayObject w, z, uh, vh, eta, xmom, *ymom;
917	int N;
918	double g, eps;
919
920	if (!PyArg_ParseTuple(args, "ddOOOOOOO",
921	&g, &eps, &w, &z, &uh, &vh, &eta,
922	&xmom, &ymom))
923	return NULL;
924
925	N = w -> dimensions[0];
926	_manning_friction_explicit(g, eps, N,
927	(double*) w -> data,
928	(double*) z -> data,
929	(double*) uh -> data,
930	(double*) vh -> data,
931	(double*) eta -> data,
932	(double*) xmom -> data,
933	(double*) ymom -> data);
934
935	return Py_BuildValue("");
936	}
937	*/
938
939
940
941	// Computational routine
942	int _extrapolate_second_order_sw(int number_of_elements,
943	double epsilon,
944	double minimum_allowed_height,
945	double beta_w,
946	double beta_w_dry,
947	double beta_uh,
948	double beta_uh_dry,
949	double beta_vh,
950	double beta_vh_dry,
951	long* surrogate_neighbours,
952	long* number_of_boundaries,
953	double* centroid_coordinates,
954	double* stage_centroid_values,
955	double* xmom_centroid_values,
956	double* ymom_centroid_values,
957	double* elevation_centroid_values,
958	double* vertex_coordinates,
959	double* stage_vertex_values,
960	double* xmom_vertex_values,
961	double* ymom_vertex_values,
962	double* elevation_vertex_values,
963	int optimise_dry_cells) {
964
965
966
967	// Local variables
968	double a, b; // Gradient vector used to calculate vertex values from centroids
969	int k,k0,k1,k2,k3,k6,coord_index,i;
970	double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle
971	double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
972	double dqv[3], qmin, qmax, hmin, hmax;
973	double hc, h0, h1, h2, beta_tmp;
974
975
976	for (k=0; k<number_of_elements; k++) {
977	k3=k*3;
978	k6=k*6;
979
980
981	if (number_of_boundaries[k]==3){
982	// No neighbours, set gradient on the triangle to zero
983
984	stage_vertex_values[k3] = stage_centroid_values[k];
985	stage_vertex_values[k3+1] = stage_centroid_values[k];
986	stage_vertex_values[k3+2] = stage_centroid_values[k];
987	xmom_vertex_values[k3] = xmom_centroid_values[k];
988	xmom_vertex_values[k3+1] = xmom_centroid_values[k];
989	xmom_vertex_values[k3+2] = xmom_centroid_values[k];
990	ymom_vertex_values[k3] = ymom_centroid_values[k];
991	ymom_vertex_values[k3+1] = ymom_centroid_values[k];
992	ymom_vertex_values[k3+2] = ymom_centroid_values[k];
993
994	continue;
995	}
996	else {
997	// Triangle k has one or more neighbours.
998	// Get centroid and vertex coordinates of the triangle
999
1000	// Get the vertex coordinates
1001	xv0 = vertex_coordinates[k6]; yv0=vertex_coordinates[k6+1];
1002	xv1 = vertex_coordinates[k6+2]; yv1=vertex_coordinates[k6+3];
1003	xv2 = vertex_coordinates[k6+4]; yv2=vertex_coordinates[k6+5];
1004
1005	// Get the centroid coordinates
1006	coord_index=2*k;
1007	x=centroid_coordinates[coord_index];
1008	y=centroid_coordinates[coord_index+1];
1009
1010	// Store x- and y- differentials for the vertices of
1011	// triangle k relative to the centroid
1012	dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
1013	dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
1014	}
1015
1016
1017
1018	if (number_of_boundaries[k]<=1){
1019
1020	//==============================================
1021	// Number of boundaries <= 1
1022	//==============================================
1023
1024
1025	// If no boundaries, auxiliary triangle is formed
1026	// from the centroids of the three neighbours
1027	// If one boundary, auxiliary triangle is formed
1028	// from this centroid and its two neighbours
1029
1030	k0=surrogate_neighbours[k3];
1031	k1=surrogate_neighbours[k3+1];
1032	k2=surrogate_neighbours[k3+2];
1033
1034	// Get the auxiliary triangle's vertex coordinates
1035	// (really the centroids of neighbouring triangles)
1036	coord_index=2*k0;
1037	x0=centroid_coordinates[coord_index];
1038	y0=centroid_coordinates[coord_index+1];
1039
1040	coord_index=2*k1;
1041	x1=centroid_coordinates[coord_index];
1042	y1=centroid_coordinates[coord_index+1];
1043
1044	coord_index=2*k2;
1045	x2=centroid_coordinates[coord_index];
1046	y2=centroid_coordinates[coord_index+1];
1047
1048	// Store x- and y- differentials for the vertices
1049	// of the auxiliary triangle
1050	dx1=x1-x0; dx2=x2-x0;
1051	dy1=y1-y0; dy2=y2-y0;
1052
1053	// Calculate 2*area of the auxiliary triangle
1054	// The triangle is guaranteed to be counter-clockwise
1055	area2 = dy2dx1 - dy1dx2;
1056
1057	// If the mesh is 'weird' near the boundary,
1058	// the triangle might be flat or clockwise:
1059	if (area2<=0) {
1060	PyErr_SetString(PyExc_RuntimeError,
1061	"shallow_water_ext.c: negative triangle area encountered");
1062	return -1;
1063	}
1064
1065	// Calculate heights of neighbouring cells
1066	hc = stage_centroid_values[k] - elevation_centroid_values[k];
1067	h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
1068	h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
1069	h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
1070	hmin = min(h0,min(h1,h2));
1071
1072
1073	if (optimise_dry_cells) {
1074	// Check if linear reconstruction is necessary for triangle k
1075	// This check will exclude dry cells.
1076
1077	hmax = max(h0,max(h1,h2));
1078	if (hmax < epsilon) {
1079	continue;
1080	}
1081	}
1082
1083
1084	//-----------------------------------
1085	// stage
1086	//-----------------------------------
1087
1088	// Calculate the difference between vertex 0 of the auxiliary
1089	// triangle and the centroid of triangle k
1090	dq0=stage_centroid_values[k0]-stage_centroid_values[k];
1091
1092	// Calculate differentials between the vertices
1093	// of the auxiliary triangle (centroids of neighbouring triangles)
1094	dq1=stage_centroid_values[k1]-stage_centroid_values[k0];
1095	dq2=stage_centroid_values[k2]-stage_centroid_values[k0];
1096
1097	// Calculate the gradient of stage on the auxiliary triangle
1098	a = dy2dq1 - dy1dq2;
1099	a /= area2;
1100	b = dx1dq2 - dx2dq1;
1101	b /= area2;
1102
1103	// Calculate provisional jumps in stage from the centroid
1104	// of triangle k to its vertices, to be limited
1105	dqv[0]=adxv0+bdyv0;
1106	dqv[1]=adxv1+bdyv1;
1107	dqv[2]=adxv2+bdyv2;
1108
1109	// Now we want to find min and max of the centroid and the
1110	// vertices of the auxiliary triangle and compute jumps
1111	// from the centroid to the min and max
1112	find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1113
1114	// Playing with dry wet interface
1115	hmin = qmin;
1116	beta_tmp = beta_w;
1117	if (hmin<minimum_allowed_height)
1118	beta_tmp = beta_w_dry;
1119
1120	// Limit the gradient
1121	limit_gradient(dqv,qmin,qmax,beta_tmp);
1122
1123	for (i=0;i<3;i++)
1124	stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
1125
1126
1127	//-----------------------------------
1128	// xmomentum
1129	//-----------------------------------
1130
1131	// Calculate the difference between vertex 0 of the auxiliary
1132	// triangle and the centroid of triangle k
1133	dq0=xmom_centroid_values[k0]-xmom_centroid_values[k];
1134
1135	// Calculate differentials between the vertices
1136	// of the auxiliary triangle
1137	dq1=xmom_centroid_values[k1]-xmom_centroid_values[k0];
1138	dq2=xmom_centroid_values[k2]-xmom_centroid_values[k0];
1139
1140	// Calculate the gradient of xmom on the auxiliary triangle
1141	a = dy2dq1 - dy1dq2;
1142	a /= area2;
1143	b = dx1dq2 - dx2dq1;
1144	b /= area2;
1145
1146	// Calculate provisional jumps in stage from the centroid
1147	// of triangle k to its vertices, to be limited
1148	dqv[0]=adxv0+bdyv0;
1149	dqv[1]=adxv1+bdyv1;
1150	dqv[2]=adxv2+bdyv2;
1151
1152	// Now we want to find min and max of the centroid and the
1153	// vertices of the auxiliary triangle and compute jumps
1154	// from the centroid to the min and max
1155	find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1156	beta_tmp = beta_uh;
1157	if (hmin<minimum_allowed_height)
1158	beta_tmp = beta_uh_dry;
1159
1160	// Limit the gradient
1161	limit_gradient(dqv,qmin,qmax,beta_tmp);
1162
1163	for (i=0;i<3;i++)
1164	xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
1165
1166
1167	//-----------------------------------
1168	// ymomentum
1169	//-----------------------------------
1170
1171	// Calculate the difference between vertex 0 of the auxiliary
1172	// triangle and the centroid of triangle k
1173	dq0=ymom_centroid_values[k0]-ymom_centroid_values[k];
1174
1175	// Calculate differentials between the vertices
1176	// of the auxiliary triangle
1177	dq1=ymom_centroid_values[k1]-ymom_centroid_values[k0];
1178	dq2=ymom_centroid_values[k2]-ymom_centroid_values[k0];
1179
1180	// Calculate the gradient of xmom on the auxiliary triangle
1181	a = dy2dq1 - dy1dq2;
1182	a /= area2;
1183	b = dx1dq2 - dx2dq1;
1184	b /= area2;
1185
1186	// Calculate provisional jumps in stage from the centroid
1187	// of triangle k to its vertices, to be limited
1188	dqv[0]=adxv0+bdyv0;
1189	dqv[1]=adxv1+bdyv1;
1190	dqv[2]=adxv2+bdyv2;
1191
1192	// Now we want to find min and max of the centroid and the
1193	// vertices of the auxiliary triangle and compute jumps
1194	// from the centroid to the min and max
1195	find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
1196	beta_tmp = beta_vh;
1197
1198	if (hmin<minimum_allowed_height)
1199	beta_tmp = beta_vh_dry;
1200
1201	// Limit the gradient
1202	limit_gradient(dqv,qmin,qmax,beta_tmp);
1203
1204	for (i=0;i<3;i++)
1205	ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
1206
1207	} // End number_of_boundaries <=1
1208	else{
1209
1210	//==============================================
1211	// Number of boundaries == 2
1212	//==============================================
1213
1214	// One internal neighbour and gradient is in direction of the neighbour's centroid
1215
1216	// Find the only internal neighbour
1217	for (k2=k3;k2<k3+3;k2++){
1218	// Find internal neighbour of triangle k
1219	// k2 indexes the edges of triangle k
1220
1221	if (surrogate_neighbours[k2]!=k)
1222	break;
1223	}
1224	if ((k2==k3+3)) {
1225	// If we didn't find an internal neighbour
1226	PyErr_SetString(PyExc_RuntimeError,
1227	"shallow_water_ext.c: Internal neighbour not found");
1228	return -1;
1229	}
1230
1231	k1=surrogate_neighbours[k2];
1232
1233	// The coordinates of the triangle are already (x,y).
1234	// Get centroid of the neighbour (x1,y1)
1235	coord_index=2*k1;
1236	x1=centroid_coordinates[coord_index];
1237	y1=centroid_coordinates[coord_index+1];
1238
1239	// Compute x- and y- distances between the centroid of
1240	// triangle k and that of its neighbour
1241	dx1=x1-x; dy1=y1-y;
1242
1243	// Set area2 as the square of the distance
1244	area2=dx1dx1+dy1dy1;
1245
1246	// Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
1247	// and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
1248	// respectively correspond to the x- and y- gradients
1249	// of the conserved quantities
1250	dx2=1.0/area2;
1251	dy2=dx2*dy1;
1252	dx2*=dx1;
1253
1254
1255
1256	//-----------------------------------
1257	// stage
1258	//-----------------------------------
1259
1260	// Compute differentials
1261	dq1=stage_centroid_values[k1]-stage_centroid_values[k];
1262
1263	// Calculate the gradient between the centroid of triangle k
1264	// and that of its neighbour
1265	a=dq1*dx2;
1266	b=dq1*dy2;
1267
1268	// Calculate provisional vertex jumps, to be limited
1269	dqv[0]=adxv0+bdyv0;
1270	dqv[1]=adxv1+bdyv1;
1271	dqv[2]=adxv2+bdyv2;
1272
1273	// Now limit the jumps
1274	if (dq1>=0.0){
1275	qmin=0.0;
1276	qmax=dq1;
1277	}
1278	else{
1279	qmin=dq1;
1280	qmax=0.0;
1281	}
1282
1283	// Limit the gradient
1284	limit_gradient(dqv,qmin,qmax,beta_w);
1285
1286	for (i=0;i<3;i++)
1287	stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
1288
1289
1290	//-----------------------------------
1291	// xmomentum
1292	//-----------------------------------
1293
1294	// Compute differentials
1295	dq1=xmom_centroid_values[k1]-xmom_centroid_values[k];
1296
1297	// Calculate the gradient between the centroid of triangle k
1298	// and that of its neighbour
1299	a=dq1*dx2;
1300	b=dq1*dy2;
1301
1302	// Calculate provisional vertex jumps, to be limited
1303	dqv[0]=adxv0+bdyv0;
1304	dqv[1]=adxv1+bdyv1;
1305	dqv[2]=adxv2+bdyv2;
1306
1307	// Now limit the jumps
1308	if (dq1>=0.0){
1309	qmin=0.0;
1310	qmax=dq1;
1311	}
1312	else{
1313	qmin=dq1;
1314	qmax=0.0;
1315	}
1316
1317	// Limit the gradient
1318	limit_gradient(dqv,qmin,qmax,beta_w);
1319
1320	for (i=0;i<3;i++)
1321	xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
1322
1323
1324	//-----------------------------------
1325	// ymomentum
1326	//-----------------------------------
1327
1328	// Compute differentials
1329	dq1=ymom_centroid_values[k1]-ymom_centroid_values[k];
1330
1331	// Calculate the gradient between the centroid of triangle k
1332	// and that of its neighbour
1333	a=dq1*dx2;
1334	b=dq1*dy2;
1335
1336	// Calculate provisional vertex jumps, to be limited
1337	dqv[0]=adxv0+bdyv0;
1338	dqv[1]=adxv1+bdyv1;
1339	dqv[2]=adxv2+bdyv2;
1340
1341	// Now limit the jumps
1342	if (dq1>=0.0){
1343	qmin=0.0;
1344	qmax=dq1;
1345	}
1346	else{
1347	qmin=dq1;
1348	qmax=0.0;
1349	}
1350
1351	// Limit the gradient
1352	limit_gradient(dqv,qmin,qmax,beta_w);
1353
1354	for (i=0;i<3;i++)
1355	ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
1356
1357	} // else [number_of_boundaries==2]
1358	} // for k=0 to number_of_elements-1
1359
1360	return 0;
1361	}
1362
1363
1364	PyObject extrapolate_second_order_sw(PyObject self, PyObject *args) {
1365	/*Compute the vertex values based on a linear reconstruction
1366	on each triangle
1367
1368	These values are calculated as follows:
1369	1) For each triangle not adjacent to a boundary, we consider the
1370	auxiliary triangle formed by the centroids of its three
1371	neighbours.
1372	2) For each conserved quantity, we integrate around the auxiliary
1373	triangle's boundary the product of the quantity and the outward
1374	normal vector. Dividing by the triangle area gives (a,b), the
1375	average of the vector (q_x,q_y) on the auxiliary triangle.
1376	We suppose that the linear reconstruction on the original
1377	triangle has gradient (a,b).
1378	3) Provisional vertex jumps dqv[0,1,2] are computed and these are
1379	then limited by calling the functions find_qmin_and_qmax and
1380	limit_gradient
1381
1382	Python call:
1383	extrapolate_second_order_sw(domain.surrogate_neighbours,
1384	domain.number_of_boundaries
1385	domain.centroid_coordinates,
1386	Stage.centroid_values
1387	Xmom.centroid_values
1388	Ymom.centroid_values
1389	domain.vertex_coordinates,
1390	Stage.vertex_values,
1391	Xmom.vertex_values,
1392	Ymom.vertex_values)
1393
1394	Post conditions:
1395	The vertices of each triangle have values from a
1396	limited linear reconstruction
1397	based on centroid values
1398
1399	*/
1400	PyArrayObject *surrogate_neighbours,
1401	*number_of_boundaries,
1402	*centroid_coordinates,
1403	*stage_centroid_values,
1404	*xmom_centroid_values,
1405	*ymom_centroid_values,
1406	*elevation_centroid_values,
1407	*vertex_coordinates,
1408	*stage_vertex_values,
1409	*xmom_vertex_values,
1410	*ymom_vertex_values,
1411	*elevation_vertex_values;
1412	PyObject domain, Tmp;
1413
1414
1415	double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry;
1416	double minimum_allowed_height, epsilon;
1417	int optimise_dry_cells, number_of_elements, e;
1418
1419	// Provisional jumps from centroids to v'tices and safety factor re limiting
1420	// by which these jumps are limited
1421	// Convert Python arguments to C
1422	if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
1423	&domain,
1424	&surrogate_neighbours,
1425	&number_of_boundaries,
1426	&centroid_coordinates,
1427	&stage_centroid_values,
1428	&xmom_centroid_values,
1429	&ymom_centroid_values,
1430	&elevation_centroid_values,
1431	&vertex_coordinates,
1432	&stage_vertex_values,
1433	&xmom_vertex_values,
1434	&ymom_vertex_values,
1435	&elevation_vertex_values,
1436	&optimise_dry_cells)) {
1437
1438	PyErr_SetString(PyExc_RuntimeError,
1439	"Input arguments to extrapolate_second_order_sw failed");
1440	return NULL;
1441	}
1442
1443	// Get the safety factor beta_w, set in the config.py file.
1444	// This is used in the limiting process
1445
1446
1447	Tmp = PyObject_GetAttrString(domain, "beta_w");
1448	if (!Tmp) {
1449	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w from domain");
1450	return NULL;
1451	}
1452	beta_w = PyFloat_AsDouble(Tmp);
1453	Py_DECREF(Tmp);
1454
1455	Tmp = PyObject_GetAttrString(domain, "beta_w_dry");
1456	if (!Tmp) {
1457	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_w_dry from domain");
1458	return NULL;
1459	}
1460	beta_w_dry = PyFloat_AsDouble(Tmp);
1461	Py_DECREF(Tmp);
1462
1463	Tmp = PyObject_GetAttrString(domain, "beta_uh");
1464	if (!Tmp) {
1465	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh from domain");
1466	return NULL;
1467	}
1468	beta_uh = PyFloat_AsDouble(Tmp);
1469	Py_DECREF(Tmp);
1470
1471	Tmp = PyObject_GetAttrString(domain, "beta_uh_dry");
1472	if (!Tmp) {
1473	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_uh_dry from domain");
1474	return NULL;
1475	}
1476	beta_uh_dry = PyFloat_AsDouble(Tmp);
1477	Py_DECREF(Tmp);
1478
1479	Tmp = PyObject_GetAttrString(domain, "beta_vh");
1480	if (!Tmp) {
1481	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh from domain");
1482	return NULL;
1483	}
1484	beta_vh = PyFloat_AsDouble(Tmp);
1485	Py_DECREF(Tmp);
1486
1487	Tmp = PyObject_GetAttrString(domain, "beta_vh_dry");
1488	if (!Tmp) {
1489	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object beta_vh_dry from domain");
1490	return NULL;
1491	}
1492	beta_vh_dry = PyFloat_AsDouble(Tmp);
1493	Py_DECREF(Tmp);
1494
1495	Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
1496	if (!Tmp) {
1497	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_height");
1498	return NULL;
1499	}
1500	minimum_allowed_height = PyFloat_AsDouble(Tmp);
1501	Py_DECREF(Tmp);
1502
1503	Tmp = PyObject_GetAttrString(domain, "epsilon");
1504	if (!Tmp) {
1505	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object epsilon");
1506	return NULL;
1507	}
1508	epsilon = PyFloat_AsDouble(Tmp);
1509	Py_DECREF(Tmp);
1510
1511
1512	// Call underlying computational routine
1513	number_of_elements = stage_centroid_values -> dimensions[0];
1514	e = _extrapolate_second_order_sw(number_of_elements,
1515	epsilon,
1516	minimum_allowed_height,
1517	beta_w,
1518	beta_w_dry,
1519	beta_uh,
1520	beta_uh_dry,
1521	beta_vh,
1522	beta_vh_dry,
1523	(long*) surrogate_neighbours -> data,
1524	(long*) number_of_boundaries -> data,
1525	(double*) centroid_coordinates -> data,
1526	(double*) stage_centroid_values -> data,
1527	(double*) xmom_centroid_values -> data,
1528	(double*) ymom_centroid_values -> data,
1529	(double*) elevation_centroid_values -> data,
1530	(double*) vertex_coordinates -> data,
1531	(double*) stage_vertex_values -> data,
1532	(double*) xmom_vertex_values -> data,
1533	(double*) ymom_vertex_values -> data,
1534	(double*) elevation_vertex_values -> data,
1535	optimise_dry_cells);
1536	if (e == -1) {
1537	// Use error string set inside computational routine
1538	return NULL;
1539	}
1540
1541
1542	return Py_BuildValue("");
1543
1544	}// extrapolate_second-order_sw
1545
1546
1547
1548
1549	PyObject rotate(PyObject self, PyObject args, PyObject kwargs) {
1550	//
1551	// r = rotate(q, normal, direction=1)
1552	//
1553	// Where q is assumed to be a Float numeric array of length 3 and
1554	// normal a Float numeric array of length 2.
1555
1556	// FIXME(Ole): I don't think this is used anymore
1557
1558	PyObject Q, Normal;
1559	PyArrayObject q, r, *normal;
1560
1561	static char *argnames[] = {"q", "normal", "direction", NULL};
1562	int dimensions[1], i, direction=1;
1563	double n1, n2;
1564
1565	// Convert Python arguments to C
1566	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO\|i", argnames,
1567	&Q, &Normal, &direction)) {
1568	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
1569	return NULL;
1570	}
1571
1572	// Input checks (convert sequences into numeric arrays)
1573	q = (PyArrayObject *)
1574	PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
1575	normal = (PyArrayObject *)
1576	PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
1577
1578
1579	if (normal -> dimensions[0] != 2) {
1580	PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
1581	return NULL;
1582	}
1583
1584	// Allocate space for return vector r (don't DECREF)
1585	dimensions[0] = 3;
1586	r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
1587
1588	// Copy
1589	for (i=0; i<3; i++) {
1590	((double ) (r -> data))[i] = ((double ) (q -> data))[i];
1591	}
1592
1593	// Get normal and direction
1594	n1 = ((double *) normal -> data)[0];
1595	n2 = ((double *) normal -> data)[1];
1596	if (direction == -1) n2 = -n2;
1597
1598	// Rotate
1599	_rotate((double *) r -> data, n1, n2);
1600
1601	// Release numeric arrays
1602	Py_DECREF(q);
1603	Py_DECREF(normal);
1604
1605	// Return result using PyArray to avoid memory leak
1606	return PyArray_Return(r);
1607	}
1608
1609
1610	// Computational function for flux computation
1611	double _compute_fluxes_central(int number_of_elements,
1612	double timestep,
1613	double epsilon,
1614	double H0,
1615	double g,
1616	long* neighbours,
1617	long* neighbour_edges,
1618	double* normals,
1619	double* edgelengths,
1620	double* radii,
1621	double* areas,
1622	long* tri_full_flag,
1623	double* stage_edge_values,
1624	double* xmom_edge_values,
1625	double* ymom_edge_values,
1626	double* bed_edge_values,
1627	double* stage_boundary_values,
1628	double* xmom_boundary_values,
1629	double* ymom_boundary_values,
1630	double* stage_explicit_update,
1631	double* xmom_explicit_update,
1632	double* ymom_explicit_update,
1633	long* already_computed_flux,
1634	double* max_speed_array,
1635	int optimise_dry_cells) {
1636
1637	// Local variables
1638	double max_speed, length, area, zl, zr;
1639	int k, i, m, n;
1640	int ki, nm=0, ki2; // Index shorthands
1641
1642	// Workspace (making them static actually made function slightly slower (Ole))
1643	double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
1644
1645	static long call=1; // Static local variable flagging already computed flux
1646
1647
1648	// Start computation
1649	call++; // Flag 'id' of flux calculation for this timestep
1650
1651	// Set explicit_update to zero for all conserved_quantities.
1652	// This assumes compute_fluxes called before forcing terms
1653	for (k=0; k<number_of_elements; k++) {
1654	stage_explicit_update[k]=0.0;
1655	xmom_explicit_update[k]=0.0;
1656	ymom_explicit_update[k]=0.0;
1657	}
1658
1659	// For all triangles
1660	for (k=0; k<number_of_elements; k++) {
1661
1662	// Loop through neighbours and compute edge flux for each
1663	for (i=0; i<3; i++) {
1664	ki = k*3+i; // Linear index (triangle k, edge i)
1665
1666	if (already_computed_flux[ki] == call)
1667	// We've already computed the flux across this edge
1668	continue;
1669
1670
1671	ql[0] = stage_edge_values[ki];
1672	ql[1] = xmom_edge_values[ki];
1673	ql[2] = ymom_edge_values[ki];
1674	zl = bed_edge_values[ki];
1675
1676	// Quantities at neighbour on nearest face
1677	n = neighbours[ki];
1678	if (n < 0) {
1679	m = -n-1; // Convert negative flag to index
1680
1681	qr[0] = stage_boundary_values[m];
1682	qr[1] = xmom_boundary_values[m];
1683	qr[2] = ymom_boundary_values[m];
1684	zr = zl; // Extend bed elevation to boundary
1685	} else {
1686	m = neighbour_edges[ki];
1687	nm = n*3+m; // Linear index (triangle n, edge m)
1688
1689	qr[0] = stage_edge_values[nm];
1690	qr[1] = xmom_edge_values[nm];
1691	qr[2] = ymom_edge_values[nm];
1692	zr = bed_edge_values[nm];
1693	}
1694
1695
1696	if (optimise_dry_cells) {
1697	// Check if flux calculation is necessary across this edge
1698	// This check will exclude dry cells.
1699	// This will also optimise cases where zl != zr as
1700	// long as both are dry
1701
1702	if ( fabs(ql[0] - zl) < epsilon &&
1703	fabs(qr[0] - zr) < epsilon ) {
1704	// Cell boundary is dry
1705
1706	already_computed_flux[ki] = call; // #k Done
1707	if (n>=0)
1708	already_computed_flux[nm] = call; // #n Done
1709
1710	max_speed = 0.0;
1711	continue;
1712	}
1713	}
1714
1715
1716	// Outward pointing normal vector (domain.normals[k, 2i:2i+2])
1717	ki2 = 2ki; //k6 + i*2
1718
1719	// Edge flux computation (triangle k, edge i)
1720	_flux_function_central(ql, qr, zl, zr,
1721	normals[ki2], normals[ki2+1],
1722	epsilon, H0, g,
1723	edgeflux, &max_speed);
1724
1725
1726	// Multiply edgeflux by edgelength
1727	length = edgelengths[ki];
1728	edgeflux[0] *= length;
1729	edgeflux[1] *= length;
1730	edgeflux[2] *= length;
1731
1732
1733	// Update triangle k with flux from edge i
1734	stage_explicit_update[k] -= edgeflux[0];
1735	xmom_explicit_update[k] -= edgeflux[1];
1736	ymom_explicit_update[k] -= edgeflux[2];
1737
1738	already_computed_flux[ki] = call; // #k Done
1739
1740
1741	// Update neighbour n with same flux but reversed sign
1742	if (n>=0) {
1743	stage_explicit_update[n] += edgeflux[0];
1744	xmom_explicit_update[n] += edgeflux[1];
1745	ymom_explicit_update[n] += edgeflux[2];
1746
1747	already_computed_flux[nm] = call; // #n Done
1748	}
1749
1750
1751	// Update timestep based on edge i and possibly neighbour n
1752	if (tri_full_flag[k] == 1) {
1753	if (max_speed > epsilon) {
1754	timestep = min(timestep, radii[k]/max_speed);
1755	if (n>=0)
1756	timestep = min(timestep, radii[n]/max_speed);
1757	}
1758	}
1759
1760	} // End edge i
1761
1762
1763	// Normalise triangle k by area and store for when all conserved
1764	// quantities get updated
1765	area = areas[k];
1766	stage_explicit_update[k] /= area;
1767	xmom_explicit_update[k] /= area;
1768	ymom_explicit_update[k] /= area;
1769
1770
1771	// Keep track of maximal speeds
1772	max_speed_array[k] = max_speed;
1773
1774	} // End triangle k
1775
1776
1777
1778	return timestep;
1779	}
1780
1781
1782	PyObject compute_fluxes_ext_central(PyObject self, PyObject *args) {
1783	/*Compute all fluxes and the timestep suitable for all volumes
1784	in domain.
1785
1786	Compute total flux for each conserved quantity using "flux_function_central"
1787
1788	Fluxes across each edge are scaled by edgelengths and summed up
1789	Resulting flux is then scaled by area and stored in
1790	explicit_update for each of the three conserved quantities
1791	stage, xmomentum and ymomentum
1792
1793	The maximal allowable speed computed by the flux_function for each volume
1794	is converted to a timestep that must not be exceeded. The minimum of
1795	those is computed as the next overall timestep.
1796
1797	Python call:
1798	domain.timestep = compute_fluxes(timestep,
1799	domain.epsilon,
1800	domain.H0,
1801	domain.g,
1802	domain.neighbours,
1803	domain.neighbour_edges,
1804	domain.normals,
1805	domain.edgelengths,
1806	domain.radii,
1807	domain.areas,
1808	tri_full_flag,
1809	Stage.edge_values,
1810	Xmom.edge_values,
1811	Ymom.edge_values,
1812	Bed.edge_values,
1813	Stage.boundary_values,
1814	Xmom.boundary_values,
1815	Ymom.boundary_values,
1816	Stage.explicit_update,
1817	Xmom.explicit_update,
1818	Ymom.explicit_update,
1819	already_computed_flux,
1820	optimise_dry_cells)
1821
1822
1823	Post conditions:
1824	domain.explicit_update is reset to computed flux values
1825	domain.timestep is set to the largest step satisfying all volumes.
1826
1827
1828	*/
1829
1830
1831	PyArrayObject neighbours, neighbour_edges,
1832	normals, edgelengths, radii, areas,
1833	*tri_full_flag,
1834	*stage_edge_values,
1835	*xmom_edge_values,
1836	*ymom_edge_values,
1837	*bed_edge_values,
1838	*stage_boundary_values,
1839	*xmom_boundary_values,
1840	*ymom_boundary_values,
1841	*stage_explicit_update,
1842	*xmom_explicit_update,
1843	*ymom_explicit_update,
1844	*already_computed_flux, //Tracks whether the flux across an edge has already been computed
1845	*max_speed_array; //Keeps track of max speeds for each triangle
1846
1847
1848	double timestep, epsilon, H0, g;
1849	int optimise_dry_cells;
1850
1851	// Convert Python arguments to C
1852	if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
1853	&timestep,
1854	&epsilon,
1855	&H0,
1856	&g,
1857	&neighbours,
1858	&neighbour_edges,
1859	&normals,
1860	&edgelengths, &radii, &areas,
1861	&tri_full_flag,
1862	&stage_edge_values,
1863	&xmom_edge_values,
1864	&ymom_edge_values,
1865	&bed_edge_values,
1866	&stage_boundary_values,
1867	&xmom_boundary_values,
1868	&ymom_boundary_values,
1869	&stage_explicit_update,
1870	&xmom_explicit_update,
1871	&ymom_explicit_update,
1872	&already_computed_flux,
1873	&max_speed_array,
1874	&optimise_dry_cells)) {
1875	PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
1876	return NULL;
1877	}
1878
1879
1880	int number_of_elements = stage_edge_values -> dimensions[0];
1881
1882	// Call underlying flux computation routine and update
1883	// the explicit update arrays
1884	timestep = _compute_fluxes_central(number_of_elements,
1885	timestep,
1886	epsilon,
1887	H0,
1888	g,
1889	(long*) neighbours -> data,
1890	(long*) neighbour_edges -> data,
1891	(double*) normals -> data,
1892	(double*) edgelengths -> data,
1893	(double*) radii -> data,
1894	(double*) areas -> data,
1895	(long*) tri_full_flag -> data,
1896	(double*) stage_edge_values -> data,
1897	(double*) xmom_edge_values -> data,
1898	(double*) ymom_edge_values -> data,
1899	(double*) bed_edge_values -> data,
1900	(double*) stage_boundary_values -> data,
1901	(double*) xmom_boundary_values -> data,
1902	(double*) ymom_boundary_values -> data,
1903	(double*) stage_explicit_update -> data,
1904	(double*) xmom_explicit_update -> data,
1905	(double*) ymom_explicit_update -> data,
1906	(long*) already_computed_flux -> data,
1907	(double*) max_speed_array -> data,
1908	optimise_dry_cells);
1909
1910	// Return updated flux timestep
1911	return Py_BuildValue("d", timestep);
1912	}
1913
1914
1915
1916	PyObject compute_fluxes_ext_kinetic(PyObject self, PyObject *args) {
1917	/*Compute all fluxes and the timestep suitable for all volumes
1918	in domain.
1919
1920	Compute total flux for each conserved quantity using "flux_function_central"
1921
1922	Fluxes across each edge are scaled by edgelengths and summed up
1923	Resulting flux is then scaled by area and stored in
1924	explicit_update for each of the three conserved quantities
1925	stage, xmomentum and ymomentum
1926
1927	The maximal allowable speed computed by the flux_function for each volume
1928	is converted to a timestep that must not be exceeded. The minimum of
1929	those is computed as the next overall timestep.
1930
1931	Python call:
1932	domain.timestep = compute_fluxes(timestep,
1933	domain.epsilon,
1934	domain.H0,
1935	domain.g,
1936	domain.neighbours,
1937	domain.neighbour_edges,
1938	domain.normals,
1939	domain.edgelengths,
1940	domain.radii,
1941	domain.areas,
1942	Stage.edge_values,
1943	Xmom.edge_values,
1944	Ymom.edge_values,
1945	Bed.edge_values,
1946	Stage.boundary_values,
1947	Xmom.boundary_values,
1948	Ymom.boundary_values,
1949	Stage.explicit_update,
1950	Xmom.explicit_update,
1951	Ymom.explicit_update,
1952	already_computed_flux)
1953
1954
1955	Post conditions:
1956	domain.explicit_update is reset to computed flux values
1957	domain.timestep is set to the largest step satisfying all volumes.
1958
1959
1960	*/
1961
1962
1963	PyArrayObject neighbours, neighbour_edges,
1964	normals, edgelengths, radii, areas,
1965	*tri_full_flag,
1966	*stage_edge_values,
1967	*xmom_edge_values,
1968	*ymom_edge_values,
1969	*bed_edge_values,
1970	*stage_boundary_values,
1971	*xmom_boundary_values,
1972	*ymom_boundary_values,
1973	*stage_explicit_update,
1974	*xmom_explicit_update,
1975	*ymom_explicit_update,
1976	*already_computed_flux; // Tracks whether the flux across an edge has already been computed
1977
1978
1979	// Local variables
1980	double timestep, max_speed, epsilon, g, H0;
1981	double normal[2], ql[3], qr[3], zl, zr;
1982	double edgeflux[3]; //Work arrays for summing up fluxes
1983
1984	int number_of_elements, k, i, m, n;
1985	int ki, nm=0, ki2; //Index shorthands
1986	static long call=1;
1987
1988
1989	// Convert Python arguments to C
1990	if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
1991	&timestep,
1992	&epsilon,
1993	&H0,
1994	&g,
1995	&neighbours,
1996	&neighbour_edges,
1997	&normals,
1998	&edgelengths, &radii, &areas,
1999	&tri_full_flag,
2000	&stage_edge_values,
2001	&xmom_edge_values,
2002	&ymom_edge_values,
2003	&bed_edge_values,
2004	&stage_boundary_values,
2005	&xmom_boundary_values,
2006	&ymom_boundary_values,
2007	&stage_explicit_update,
2008	&xmom_explicit_update,
2009	&ymom_explicit_update,
2010	&already_computed_flux)) {
2011	PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
2012	return NULL;
2013	}
2014	number_of_elements = stage_edge_values -> dimensions[0];
2015	call++;//a static local variable to which already_computed_flux is compared
2016	//set explicit_update to zero for all conserved_quantities.
2017	//This assumes compute_fluxes called before forcing terms
2018	for (k=0; k<number_of_elements; k++) {
2019	((double *) stage_explicit_update -> data)[k]=0.0;
2020	((double *) xmom_explicit_update -> data)[k]=0.0;
2021	((double *) ymom_explicit_update -> data)[k]=0.0;
2022	}
2023	//Loop through neighbours and compute edge flux for each
2024	for (k=0; k<number_of_elements; k++) {
2025	for (i=0; i<3; i++) {
2026	ki = k*3+i;
2027	if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
2028	continue;
2029	ql[0] = ((double *) stage_edge_values -> data)[ki];
2030	ql[1] = ((double *) xmom_edge_values -> data)[ki];
2031	ql[2] = ((double *) ymom_edge_values -> data)[ki];
2032	zl = ((double *) bed_edge_values -> data)[ki];
2033
2034	//Quantities at neighbour on nearest face
2035	n = ((long *) neighbours -> data)[ki];
2036	if (n < 0) {
2037	m = -n-1; //Convert negative flag to index
2038	qr[0] = ((double *) stage_boundary_values -> data)[m];
2039	qr[1] = ((double *) xmom_boundary_values -> data)[m];
2040	qr[2] = ((double *) ymom_boundary_values -> data)[m];
2041	zr = zl; //Extend bed elevation to boundary
2042	} else {
2043	m = ((long *) neighbour_edges -> data)[ki];
2044	nm = n*3+m;
2045	qr[0] = ((double *) stage_edge_values -> data)[nm];
2046	qr[1] = ((double *) xmom_edge_values -> data)[nm];
2047	qr[2] = ((double *) ymom_edge_values -> data)[nm];
2048	zr = ((double *) bed_edge_values -> data)[nm];
2049	}
2050	// Outward pointing normal vector
2051	// normal = domain.normals[k, 2i:2i+2]
2052	ki2 = 2ki; //k6 + i*2
2053	normal[0] = ((double *) normals -> data)[ki2];
2054	normal[1] = ((double *) normals -> data)[ki2+1];
2055	//Edge flux computation
2056	flux_function_kinetic(ql, qr, zl, zr,
2057	normal[0], normal[1],
2058	epsilon, H0, g,
2059	edgeflux, &max_speed);
2060	//update triangle k
2061	((long *) already_computed_flux->data)[ki]=call;
2062	((double ) stage_explicit_update -> data)[k] -= edgeflux[0]((double *) edgelengths -> data)[ki];
2063	((double ) xmom_explicit_update -> data)[k] -= edgeflux[1]((double *) edgelengths -> data)[ki];
2064	((double ) ymom_explicit_update -> data)[k] -= edgeflux[2]((double *) edgelengths -> data)[ki];
2065	//update the neighbour n
2066	if (n>=0){
2067	((long *) already_computed_flux->data)[nm]=call;
2068	((double ) stage_explicit_update -> data)[n] += edgeflux[0]((double *) edgelengths -> data)[nm];
2069	((double ) xmom_explicit_update -> data)[n] += edgeflux[1]((double *) edgelengths -> data)[nm];
2070	((double ) ymom_explicit_update -> data)[n] += edgeflux[2]((double *) edgelengths -> data)[nm];
2071	}
2072	///for (j=0; j<3; j++) {
2073	///flux[j] -= edgeflux[j]((double ) edgelengths -> data)[ki];
2074	///}
2075	//Update timestep
2076	//timestep = min(timestep, domain.radii[k]/max_speed)
2077	//FIXME: SR Add parameter for CFL condition
2078	if ( ((long *) tri_full_flag -> data)[k] == 1) {
2079	if (max_speed > epsilon) {
2080	timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
2081	//maxspeed in flux_function is calculated as max(\|u+a\|,\|u-a\|)
2082	if (n>=0)
2083	timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
2084	}
2085	}
2086	} // end for i
2087	//Normalise by area and store for when all conserved
2088	//quantities get updated
2089	((double ) stage_explicit_update -> data)[k] /= ((double ) areas -> data)[k];
2090	((double ) xmom_explicit_update -> data)[k] /= ((double ) areas -> data)[k];
2091	((double ) ymom_explicit_update -> data)[k] /= ((double ) areas -> data)[k];
2092	} //end for k
2093	return Py_BuildValue("d", timestep);
2094	}
2095
2096	PyObject protect(PyObject self, PyObject *args) {
2097	//
2098	// protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
2099
2100
2101	PyArrayObject
2102	*wc, //Stage at centroids
2103	*zc, //Elevation at centroids
2104	*xmomc, //Momentums at centroids
2105	*ymomc;
2106
2107
2108	int N;
2109	double minimum_allowed_height, maximum_allowed_speed, epsilon;
2110
2111	// Convert Python arguments to C
2112	if (!PyArg_ParseTuple(args, "dddOOOO",
2113	&minimum_allowed_height,
2114	&maximum_allowed_speed,
2115	&epsilon,
2116	&wc, &zc, &xmomc, &ymomc)) {
2117	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
2118	return NULL;
2119	}
2120
2121	N = wc -> dimensions[0];
2122
2123	_protect(N,
2124	minimum_allowed_height,
2125	maximum_allowed_speed,
2126	epsilon,
2127	(double*) wc -> data,
2128	(double*) zc -> data,
2129	(double*) xmomc -> data,
2130	(double*) ymomc -> data);
2131
2132	return Py_BuildValue("");
2133	}
2134
2135
2136
2137	PyObject balance_deep_and_shallow(PyObject self, PyObject *args) {
2138	// Compute linear combination between stage as computed by
2139	// gradient-limiters limiting using w, and stage computed by
2140	// gradient-limiters limiting using h (h-limiter).
2141	// The former takes precedence when heights are large compared to the
2142	// bed slope while the latter takes precedence when heights are
2143	// relatively small. Anything in between is computed as a balanced
2144	// linear combination in order to avoid numerical disturbances which
2145	// would otherwise appear as a result of hard switching between
2146	// modes.
2147	//
2148	// balance_deep_and_shallow(beta_h, wc, zc, wv, zv,
2149	// xmomc, ymomc, xmomv, ymomv)
2150
2151
2152	PyArrayObject
2153	*wc, //Stage at centroids
2154	*zc, //Elevation at centroids
2155	*wv, //Stage at vertices
2156	*zv, //Elevation at vertices
2157	*hvbar, //h-Limited depths at vertices
2158	*xmomc, //Momentums at centroids and vertices
2159	*ymomc,
2160	*xmomv,
2161	*ymomv;
2162
2163	PyObject domain, Tmp;
2164
2165	double alpha_balance = 2.0;
2166	double H0, beta_h;
2167
2168	int N, tight_slope_limiters; //, err;
2169
2170	// Convert Python arguments to C
2171	if (!PyArg_ParseTuple(args, "OdOOOOOOOOO",
2172	&domain,
2173	&beta_h,
2174	&wc, &zc,
2175	&wv, &zv, &hvbar,
2176	&xmomc, &ymomc, &xmomv, &ymomv)) {
2177	PyErr_SetString(PyExc_RuntimeError,
2178	"shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
2179	return NULL;
2180	}
2181
2182
2183	// FIXME (Ole): I tested this without GetAttrString and got time down
2184	// marginally from 4.0s to 3.8s. Consider passing everything in
2185	// through ParseTuple and profile.
2186
2187	// Pull out parameters
2188	Tmp = PyObject_GetAttrString(domain, "alpha_balance");
2189	if (!Tmp) {
2190	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
2191	return NULL;
2192	}
2193	alpha_balance = PyFloat_AsDouble(Tmp);
2194	Py_DECREF(Tmp);
2195
2196
2197	Tmp = PyObject_GetAttrString(domain, "H0");
2198	if (!Tmp) {
2199	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
2200	return NULL;
2201	}
2202	H0 = PyFloat_AsDouble(Tmp);
2203	Py_DECREF(Tmp);
2204
2205
2206	Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
2207	if (!Tmp) {
2208	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object tight_slope_limiters from domain");
2209	return NULL;
2210	}
2211	tight_slope_limiters = PyInt_AsLong(Tmp);
2212	Py_DECREF(Tmp);
2213
2214
2215	N = wc -> dimensions[0];
2216	_balance_deep_and_shallow(N,
2217	beta_h,
2218	(double*) wc -> data,
2219	(double*) zc -> data,
2220	(double*) wv -> data,
2221	(double*) zv -> data,
2222	(double*) hvbar -> data,
2223	(double*) xmomc -> data,
2224	(double*) ymomc -> data,
2225	(double*) xmomv -> data,
2226	(double*) ymomv -> data,
2227	H0,
2228	(int) tight_slope_limiters,
2229	alpha_balance);
2230
2231
2232	return Py_BuildValue("");
2233	}
2234
2235
2236
2237	PyObject h_limiter(PyObject self, PyObject *args) {
2238
2239	PyObject domain, Tmp;
2240	PyArrayObject
2241	hv, hc, //Depth at vertices and centroids
2242	*hvbar, //Limited depth at vertices (return values)
2243	*neighbours;
2244
2245	int k, i, n, N, k3;
2246	int dimensions[2];
2247	double beta_h; // Safety factor (see config.py)
2248	double hmin, hmax, hn;
2249
2250	// Convert Python arguments to C
2251	if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
2252	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not parse input arguments");
2253	return NULL;
2254	}
2255
2256	neighbours = get_consecutive_array(domain, "neighbours");
2257
2258	//Get safety factor beta_h
2259	Tmp = PyObject_GetAttrString(domain, "beta_h");
2260	if (!Tmp) {
2261	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter could not obtain object beta_h from domain");
2262	return NULL;
2263	}
2264	beta_h = PyFloat_AsDouble(Tmp);
2265
2266	Py_DECREF(Tmp);
2267
2268	N = hc -> dimensions[0];
2269
2270	//Create hvbar
2271	dimensions[0] = N;
2272	dimensions[1] = 3;
2273	hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
2274
2275
2276	//Find min and max of this and neighbour's centroid values
2277	hmin = malloc(N * sizeof(double));
2278	hmax = malloc(N * sizeof(double));
2279	for (k=0; k<N; k++) {
2280	k3=k*3;
2281
2282	hmin[k] = ((double*) hc -> data)[k];
2283	hmax[k] = hmin[k];
2284
2285	for (i=0; i<3; i++) {
2286	n = ((long*) neighbours -> data)[k3+i];
2287
2288	// Initialise hvbar with values from hv
2289	((double) hvbar -> data)[k3+i] = ((double) hv -> data)[k3+i];
2290
2291	if (n >= 0) {
2292	hn = ((double*) hc -> data)[n]; // Neighbour's centroid value
2293
2294	hmin[k] = min(hmin[k], hn);
2295	hmax[k] = max(hmax[k], hn);
2296	}
2297	}
2298	}
2299
2300	// Call underlying standard routine
2301	_limit(N, beta_h, (double) hc -> data, (double) hvbar -> data, hmin, hmax);
2302
2303	// // //Py_DECREF(domain); //FIXME: Necessary?
2304	free(hmin);
2305	free(hmax);
2306
2307	// Return result using PyArray to avoid memory leak
2308	return PyArray_Return(hvbar);
2309	}
2310
2311	PyObject h_limiter_sw(PyObject self, PyObject *args) {
2312	//a faster version of h_limiter above
2313	PyObject domain, Tmp;
2314	PyArrayObject
2315	hv, hc, //Depth at vertices and centroids
2316	*hvbar, //Limited depth at vertices (return values)
2317	*neighbours;
2318
2319	int k, i, N, k3,k0,k1,k2;
2320	int dimensions[2];
2321	double beta_h; //Safety factor (see config.py)
2322	double hmin, hmax, dh[3];
2323	// Convert Python arguments to C
2324	if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv)) {
2325	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not parse input arguments");
2326	return NULL;
2327	}
2328	neighbours = get_consecutive_array(domain, "neighbours");
2329
2330	//Get safety factor beta_h
2331	Tmp = PyObject_GetAttrString(domain, "beta_h");
2332	if (!Tmp) {
2333	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: h_limiter_sw could not obtain object beta_h from domain");
2334	return NULL;
2335	}
2336	beta_h = PyFloat_AsDouble(Tmp);
2337
2338	Py_DECREF(Tmp);
2339
2340	N = hc -> dimensions[0];
2341
2342	//Create hvbar
2343	dimensions[0] = N;
2344	dimensions[1] = 3;
2345	hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
2346	for (k=0;k<N;k++){
2347	k3=k*3;
2348	//get the ids of the neighbours
2349	k0 = ((long*) neighbours -> data)[k3];
2350	k1 = ((long*) neighbours -> data)[k3+1];
2351	k2 = ((long*) neighbours -> data)[k3+2];
2352	//set hvbar provisionally
2353	for (i=0;i<3;i++){
2354	((double) hvbar -> data)[k3+i] = ((double) hv -> data)[k3+i];
2355	dh[i]=((double) hvbar -> data)[k3+i]-((double) hc -> data)[k];
2356	}
2357	hmin=((double*) hc -> data)[k];
2358	hmax=hmin;
2359	if (k0>=0){
2360	hmin=min(hmin,((double*) hc -> data)[k0]);
2361	hmax=max(hmax,((double*) hc -> data)[k0]);
2362	}
2363	if (k1>=0){
2364	hmin=min(hmin,((double*) hc -> data)[k1]);
2365	hmax=max(hmax,((double*) hc -> data)[k1]);
2366	}
2367	if (k2>=0){
2368	hmin=min(hmin,((double*) hc -> data)[k2]);
2369	hmax=max(hmax,((double*) hc -> data)[k2]);
2370	}
2371	hmin-=((double*) hc -> data)[k];
2372	hmax-=((double*) hc -> data)[k];
2373	limit_gradient(dh,hmin,hmax,beta_h);
2374	for (i=0;i<3;i++)
2375	((double) hvbar -> data)[k3+i] = ((double) hc -> data)[k]+dh[i];
2376	}
2377	return PyArray_Return(hvbar);
2378	}
2379
2380	PyObject assign_windfield_values(PyObject self, PyObject *args) {
2381	//
2382	// assign_windfield_values(xmom_update, ymom_update,
2383	// s_vec, phi_vec, self.const)
2384
2385
2386
2387	PyArrayObject //(one element per triangle)
2388	*s_vec, //Speeds
2389	*phi_vec, //Bearings
2390	*xmom_update, //Momentum updates
2391	*ymom_update;
2392
2393
2394	int N;
2395	double cw;
2396
2397	// Convert Python arguments to C
2398	if (!PyArg_ParseTuple(args, "OOOOd",
2399	&xmom_update,
2400	&ymom_update,
2401	&s_vec, &phi_vec,
2402	&cw)) {
2403	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
2404	return NULL;
2405	}
2406
2407
2408	N = xmom_update -> dimensions[0];
2409
2410	_assign_wind_field_values(N,
2411	(double*) xmom_update -> data,
2412	(double*) ymom_update -> data,
2413	(double*) s_vec -> data,
2414	(double*) phi_vec -> data,
2415	cw);
2416
2417	return Py_BuildValue("");
2418	}
2419
2420
2421
2422
2423	//-------------------------------
2424	// Method table for python module
2425	//-------------------------------
2426	static struct PyMethodDef MethodTable[] = {
2427	/* The cast of the function is necessary since PyCFunction values
2428	* only take two PyObject* parameters, and rotate() takes
2429	* three.
2430	*/
2431
2432	{"rotate", (PyCFunction)rotate, METH_VARARGS \| METH_KEYWORDS, "Print out"},
2433	{"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
2434	{"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
2435	{"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
2436	{"gravity", gravity, METH_VARARGS, "Print out"},
2437	{"manning_friction", manning_friction, METH_VARARGS, "Print out"},
2438	{"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},
2439	{"balance_deep_and_shallow", balance_deep_and_shallow,
2440	METH_VARARGS, "Print out"},
2441	{"h_limiter", h_limiter,
2442	METH_VARARGS, "Print out"},
2443	{"h_limiter_sw", h_limiter_sw,
2444	METH_VARARGS, "Print out"},
2445	{"protect", protect, METH_VARARGS \| METH_KEYWORDS, "Print out"},
2446	{"assign_windfield_values", assign_windfield_values,
2447	METH_VARARGS \| METH_KEYWORDS, "Print out"},
2448	{NULL, NULL}
2449	};
2450
2451	// Module initialisation
2452	void initshallow_water_ext(void){
2453	Py_InitModule("shallow_water_ext", MethodTable);
2454
2455	import_array(); // Necessary for handling of NumPY structures
2456	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: