Context Navigation

source: trunk/anuga_work/development/gareth/balanced_basic/swb2_domain_ext.c @ 8302

Last change on this file since 8302 was 8302, checked in by davies, 12 years ago
Adding files for swb2
File size: 39.4 KB

Line
1	// Python - C extension module for shallow_water.py
2	//
3	// To compile (Python2.6):
4	// gcc -c swb2_domain_ext.c -I/usr/include/python2.6 -o domain_ext.o -Wall -O
5	// gcc -shared swb2_domain_ext.o -o swb2_domain_ext.so
6	//
7	// or use python compile.py
8	//
9	// See the module swb_domain.py for more documentation on
10	// how to use this module
11	//
12	//
13	// Stephen Roberts, ANU 2009
14	// Ole Nielsen, GA 2004
15	// Gareth Davies, GA 2011
16
17	#include "Python.h"
18	#include "numpy/arrayobject.h"
19	#include "math.h"
20	#include <stdio.h>
21	//#include "numpy_shim.h"
22
23	// Shared code snippets
24	#include "util_ext.h"
25
26
27	const double pi = 3.14159265358979;
28
29	// Computational function for rotation
30	int _rotate(double *q, double n1, double n2) {
31	/*Rotate the last 2 coordinates of q (q[1], q[2])
32	from x,y coordinates to coordinates based on normal vector (n1, n2).
33
34	Result is returned in array 2x1 r
35	To rotate in opposite direction, call rotate with (q, n1, -n2)
36
37	Contents of q are changed by this function */
38
39
40	double q1, q2;
41
42	// Shorthands
43	q1 = q[1]; // x coordinate
44	q2 = q[2]; // y coordinate
45
46	// Rotate
47	q[1] = n1q1 + n2q2;
48	q[2] = -n2q1 + n1q2;
49
50	return 0;
51	}
52
53	// Function to obtain speed from momentum and depth.
54	// This is used by flux functions
55	// Input parameters uh and h may be modified by this function.
56	// Tried to inline, but no speedup was achieved 27th May 2009 (Ole)
57	//static inline double _compute_speed(double *uh,
58	double _compute_speed(double *uh,
59	double *h,
60	double epsilon,
61	double h0,
62	double limiting_threshold) {
63
64	double u;
65
66	if (*h < limiting_threshold) {
67	// Apply limiting of speeds according to the ANUGA manual
68	if (*h < epsilon) {
69	*h = 0.0; // Could have been negative
70	u = 0.0;
71	} else {
72	u = uh/(h + h0/ *h);
73	}
74
75
76	// Adjust momentum to be consistent with speed
77	uh = u *h;
78	} else {
79	// We are in deep water - no need for limiting
80	u = uh/ h;
81	}
82
83	return u;
84	}
85
86	// Innermost flux function (using stage w=z+h)
87	int _flux_function_central(double q_left, double q_right,
88	double z_left, double z_right,
89	double n1, double n2,
90	double epsilon,
91	double h0,
92	double limiting_threshold,
93	double g,
94	double edgeflux, double max_speed,
95	double *pressure_flux)
96	{
97
98	/*Compute fluxes between volumes for the shallow water wave equation
99	cast in terms of the 'stage', w = h+z using
100	the 'central scheme' as described in
101
102	Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
103	Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
104	Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
105
106	The implemented formula is given in equation (3.15) on page 714
107	*/
108
109	int i;
110
111	double w_left, h_left, uh_left, vh_left, u_left;
112	double w_right, h_right, uh_right, vh_right, u_right;
113	double s_min, s_max, soundspeed_left, soundspeed_right;
114	double denom, inverse_denominator, z;
115	double uint, t1, t2, t3;
116	// Workspace (allocate once, use many)
117	static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
118
119	// Copy conserved quantities to protect from modification
120	q_left_rotated[0] = q_left[0];
121	q_right_rotated[0] = q_right[0];
122	q_left_rotated[1] = q_left[1];
123	q_right_rotated[1] = q_right[1];
124	q_left_rotated[2] = q_left[2];
125	q_right_rotated[2] = q_right[2];
126
127	// Align x- and y-momentum with x-axis
128	_rotate(q_left_rotated, n1, n2);
129	_rotate(q_right_rotated, n1, n2);
130
131	z = 0.5*(z_left + z_right); // Average elevation values.
132	// Even though this will nominally allow
133	// for discontinuities in the elevation data,
134	// there is currently no numerical support for
135	// this so results may be strange near
136	// jumps in the bed.
137
138	// Compute speeds in x-direction
139	w_left = q_left_rotated[0];
140	h_left = w_left - z;
141	uh_left = q_left_rotated[1];
142	u_left = _compute_speed(&uh_left, &h_left,
143	epsilon, h0, limiting_threshold);
144
145	w_right = q_right_rotated[0];
146	h_right = w_right - z;
147	uh_right = q_right_rotated[1];
148	u_right = _compute_speed(&uh_right, &h_right,
149	epsilon, h0, limiting_threshold);
150
151	// Momentum in y-direction
152	vh_left = q_left_rotated[2];
153	vh_right = q_right_rotated[2];
154
155	// Limit y-momentum if necessary
156	// Leaving this out, improves speed significantly (Ole 27/5/2009)
157	// All validation tests pass, so do we really need it anymore?
158	_compute_speed(&vh_left, &h_left,
159	epsilon, h0, limiting_threshold);
160	_compute_speed(&vh_right, &h_right,
161	epsilon, h0, limiting_threshold);
162
163	// Maximal and minimal wave speeds
164	soundspeed_left = sqrt(g*h_left);
165	soundspeed_right = sqrt(g*h_right);
166
167	// Code to use fast square root optimisation if desired.
168	// Timings on AMD 64 for the Okushiri profile gave the following timings
169	//
170	// SQRT Total Flux
171	//=============================
172	//
173	// Ref 405s 152s
174	// Fast (dbl) 453s 173s
175	// Fast (sng) 437s 171s
176	//
177	// Consequently, there is currently (14/5/2009) no reason to use this
178	// approximation.
179
180	//soundspeed_left = fast_squareroot_approximation(g*h_left);
181	//soundspeed_right = fast_squareroot_approximation(g*h_right);
182
183	s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
184	if (s_max < 0.0)
185	{
186	s_max = 0.0;
187	}
188
189	s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
190	if (s_min > 0.0)
191	{
192	s_min = 0.0;
193	}
194
195	// Flux formulas
196	flux_left[0] = u_left*h_left;
197	flux_left[1] = u_leftuh_left + 0.5gh_lefth_left;
198	flux_left[2] = u_left*vh_left;
199
200	flux_right[0] = u_right*h_right;
201	flux_right[1] = u_rightuh_right + 0.5gh_righth_right;
202	flux_right[2] = u_right*vh_right;
203
204	// Flux computation
205	denom = s_max - s_min;
206	if (denom < epsilon)
207	{ // FIXME (Ole): Try using h0 here
208	memset(edgeflux, 0, 3*sizeof(double));
209	//for(i=0;i<3;i++){
210	// edgeflux[i] = _minmod(flux_left[i], flux_right[i]);
211	//}
212	*max_speed = 0.0;
213	}
214	else
215	{
216	inverse_denominator = 1.0/denom;
217	for (i = 0; i < 3; i++)
218	{
219	// Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
220	// Physics 2:141-163
221	//uint = (s_maxq_right_rotated[i] - s_minq_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
222	//t1 = (q_right_rotated[i] - uint);
223	//t2 = (-q_left_rotated[i] + uint);
224	//t3 = _minmod(t1,t2);
225	t3 = 0.0;
226	edgeflux[i] = s_maxflux_left[i] - s_minflux_right[i];
227	edgeflux[i] += s_maxs_min(q_right_rotated[i] - q_left_rotated[i] - t3);
228	//t1 = s_maxs_min(q_right_rotated[i] - q_left_rotated[i] - t3);
229	//if(fabs(edgeflux[i])>fabs(t1)){
230	// edgeflux[i]+=t1;
231	//}else{
232	// edgeflux[i]+=fabs(edgeflux[i])*copysign(1.0,t1) ;
233	//}
234
235	edgeflux[i] *= inverse_denominator;
236	}
237
238	pressure_flux = 0.5g( s_maxh_lefth_left -s_minh_righth_right)inverse_denominator;
239	//if(edgeflux[0]<0.0){
240	// //edgeflux[2] = edgeflux[0]0.5(vh_right/h_right + vh_left/h_left);
241	// edgeflux[2] = edgeflux[0]*vh_right/h_right;
242	//}else{
243	// edgeflux[2] = edgeflux[0]*vh_left/h_left;
244	//}
245
246	// Maximal wavespeed
247	*max_speed = max(fabs(s_max), fabs(s_min));
248
249	// Rotate back
250	_rotate(edgeflux, n1, -n2);
251	}
252
253	return 0;
254	}
255
256
257	// Computational function for flux computation
258	double _compute_fluxes_central(int number_of_elements,
259	double timestep,
260	double epsilon,
261	double H0,
262	double g,
263	long* neighbours,
264	long* neighbour_edges,
265	double* normals,
266	double* edgelengths,
267	double* radii,
268	double* areas,
269	long* tri_full_flag,
270	double* stage_edge_values,
271	double* xmom_edge_values,
272	double* ymom_edge_values,
273	double* bed_edge_values,
274	double* stage_boundary_values,
275	double* xmom_boundary_values,
276	double* ymom_boundary_values,
277	double* stage_explicit_update,
278	double* xmom_explicit_update,
279	double* ymom_explicit_update,
280	long* already_computed_flux,
281	double* max_speed_array,
282	int optimise_dry_cells,
283	double* stage_centroid_values,
284	double* bed_centroid_values,
285	double* bed_vertex_values) {
286	// Local variables
287	double max_speed, length, inv_area, zl, zr;
288	double h0 = H0*H0; // This ensures a good balance when h approaches H0.
289
290	double limiting_threshold = 10 * H0; // Avoid applying limiter below this
291	// threshold for performance reasons.
292	// See ANUGA manual under flux limiting
293	int k, i, m, n,j;
294	int ki, nm = 0, ki2; // Index shorthands
295
296	// Workspace (making them static actually made function slightly slower (Ole))
297	double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
298	double stage_edges[3];//Work array
299	double bedslope_work;
300	int neighbours_wet[3];//Work array
301	int useint;
302	double usedouble, scale_factor_shallow, bedtop, bedbot, pressure_flux;
303	//double depthtemp, max_bed_edge_value;
304	static long call = 1; // Static local variable flagging already computed flux
305
306	// Start computation
307	call++; // Flag 'id' of flux calculation for this timestep
308
309	// Set explicit_update to zero for all conserved_quantities.
310	// This assumes compute_fluxes called before forcing terms
311	memset((char) stage_explicit_update, 0, number_of_elements sizeof (double));
312	memset((char) xmom_explicit_update, 0, number_of_elements sizeof (double));
313	memset((char) ymom_explicit_update, 0, number_of_elements sizeof (double));
314
315	// For all triangles
316	for (k = 0; k < number_of_elements; k++) {
317	// Loop through neighbours and compute edge flux for each
318	for (i = 0; i < 3; i++) {
319	ki = k * 3 + i; // Linear index to edge i of triangle k
320
321	if (already_computed_flux[ki] == call) {
322	// We've already computed the flux across this edge
323	continue;
324	}
325
326	// Get left hand side values from triangle k, edge i
327	ql[0] = stage_edge_values[ki];
328	ql[1] = xmom_edge_values[ki];
329	ql[2] = ymom_edge_values[ki];
330	zl = bed_edge_values[ki];
331
332	// Get right hand side values either from neighbouring triangle
333	// or from boundary array (Quantities at neighbour on nearest face).
334	n = neighbours[ki];
335	if (n < 0) {
336	// Neighbour is a boundary condition
337	m = -n - 1; // Convert negative flag to boundary index
338
339	qr[0] = stage_boundary_values[m];
340	qr[1] = xmom_boundary_values[m];
341	qr[2] = ymom_boundary_values[m];
342	zr = zl; // Extend bed elevation to boundary
343	}
344	else {
345	// Neighbour is a real triangle
346	m = neighbour_edges[ki];
347	nm = n * 3 + m; // Linear index (triangle n, edge m)
348
349	qr[0] = stage_edge_values[nm];
350	qr[1] = xmom_edge_values[nm];
351	qr[2] = ymom_edge_values[nm];
352	zr = bed_edge_values[nm];
353	}
354
355	// Now we have values for this edge - both from left and right side.
356
357	if (optimise_dry_cells) {
358	// Check if flux calculation is necessary across this edge
359	// This check will exclude dry cells.
360	// This will also optimise cases where zl != zr as
361	// long as both are dry
362
363	if ((ql[0] - zl) < epsilon &&
364	(qr[0] - zr) < epsilon) {
365	// Cell boundary is dry
366
367	already_computed_flux[ki] = call; // #k Done
368	if (n >= 0) {
369	already_computed_flux[nm] = call; // #n Done
370	}
371
372	max_speed = 0.0;
373	continue;
374	}
375	}
376
377
378	if (fabs(zl-zr)>1.0e-10) {
379	report_python_error(AT,"Discontinuous Elevation");
380	return 0.0;
381	}
382
383	// If both centroids are dry, then the cells should not exchange flux
384	// NOTE: This also means no bedslope term -- which is arguably
385	// appropriate, since the depth is zero at the cell centroid
386	if( (stage_centroid_values[k]-bed_centroid_values[k]<=0.0)){
387	if(n>0){
388	if(stage_centroid_values[n]-bed_centroid_values[n]<=0.0){
389	continue;
390	}
391	}else{
392	// Treat boundaries -- if the boundary stage < zl, there is
393	// no flux
394	if(qr[0]< zl){
395	continue;
396	}
397	}
398	}
399
400	//If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
401	// unless the local centroid value is smaller
402	if(n>=0){
403	if((stage_centroid_values[k]-bed_centroid_values[k]<= 0.0)){
404	ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
405	}
406	if(stage_centroid_values[n]-bed_centroid_values[n]<= 0.0){
407	qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
408	}
409	}else{
410	// Treat the boundary case
411	if((stage_centroid_values[k]-bed_centroid_values[k]<=0.0)){
412	ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
413	}
414
415	}
416
417	// Outward pointing normal vector (domain.normals[k, 2i:2i+2])
418	ki2 = 2 * ki; //k6 + i2
419
420	// Edge flux computation (triangle k, edge i)
421	_flux_function_central(ql, qr, zl, zr,
422	normals[ki2], normals[ki2 + 1],
423	epsilon, h0, limiting_threshold, g,
424	edgeflux, &max_speed, &pressure_flux);
425	//_flux_function_fraccaro(ql, qr, zl, zr,
426	// normals[ki2], normals[ki2 + 1],
427	// epsilon, h0, limiting_threshold, g,
428	// edgeflux, &max_speed);
429
430	//Find the bed vertex values associated with the edge
431	//if(i==0){
432	// bedtop = max(bed_vertex_values[3k+1], bed_vertex_values[3k+2]);
433	// bedbot = min(bed_vertex_values[3k+1], bed_vertex_values[3k+2]);
434	//}else{
435	// if(i==1){
436	// bedtop = max(bed_vertex_values[3k+0], bed_vertex_values[3k+2]);
437	// bedbot = min(bed_vertex_values[3k+0], bed_vertex_values[3k+2]);
438	// }else{
439	// // i==2
440	// bedtop = max(bed_vertex_values[3k+1], bed_vertex_values[3k+0]);
441	// bedbot = min(bed_vertex_values[3k+1], bed_vertex_values[3k+0]);
442	// }
443	//
444	//}
445
446	//if(fabs(bedtop-bedbot)<1.0e-10){
447	// scale_factor_shallow = 1.0;
448	//}else{
449	// scale_factor_shallow = min( (bedtop - bedbot)/(min(stage_centroid_values[k], bedtop) - bedbot),10.0);
450	//}
451
452	// Multiply edgeflux by edgelength
453	length = edgelengths[ki];
454	//length = edgelengths[ki]*scale_factor_shallow;
455	edgeflux[0] *= length;
456	edgeflux[1] *= length;
457	edgeflux[2] *= length;
458
459	// To correctly implement reflective boundary conditions would be something like this?
460	//if(n<0){
461	// _rotate(edgeflux,normals[ki2],normals[ki2+1]);
462	// edgeflux[1] = 0.0;
463	// _rotate(edgeflux,-normals[ki2+1],normals[ki2]);
464	//}
465
466	// Update triangle k with flux from edge i
467	stage_explicit_update[k] -= edgeflux[0];
468	xmom_explicit_update[k] -= edgeflux[1];
469	ymom_explicit_update[k] -= edgeflux[2];
470
471	// Calculate the bed slope term, using the 'divergence form' from
472	// Alessandro Valiani and Lorenzo Begnudelli, Divergence Form for Bed Slope Source Term in Shallow Water Equations
473	// Journal of Hydraulic Engineering, 2006, 132:652-665
474	if((stage_centroid_values[k]>bed_centroid_values[k]+0.0)){
475	bedslope_work = -0.5gmax(stage_centroid_values[k]-zl,0.)(stage_centroid_values[k]-zl)length;
476	xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
477	ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
478	}else{
479	// Treat nearly dry cells
480	//bedslope_work =-0.5gmax(ql[0]-zl,0.)(ql[0]-zl)length; //
481	bedslope_work = -pressure_flux*length;
482	xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
483	ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
484	}
485	// The triangle is (partially?) dry
486	// Find the average stage of the wet neighbours edges -- set this as the within triangle stage for the bedslope term
487	// FIXME: Note that we keep recomputing this -- could improve efficiency
488	//for(j=0;j<3;j++){
489	// // For each edge, check if each neighbouring triangle is wet
490	// useint = neighbours[3*k+j];
491	// if(useint>0){
492	// // neighbour is not a boundary condition
493	// if(stage_centroid_values[useint] > bed_centroid_values[useint]){
494	// //wet
495	// neighbours_wet[j]=1;
496	// useint = useint3 + neighbour_edges[3k+j]; // Index of neighbouring edge
497	// stage_edges[j] = stage_edge_values[useint]; // stage of neighbouring edge
498	// }else{
499	// //dry
500	// neighbours_wet[j]=0;
501	// stage_edges[j]=0.0; // Arbitrary
502	// }
503	// }else{
504	// // neighbour is a boundary condition, ASSUME IT IS WET
505	// useint = -useint -1;
506	// neighbours_wet[j]=1;
507	// stage_edges[j] = stage_boundary_values[useint];
508	// }
509	//}
510	////calculate average neighbouring edge stage, only considering wet triangles
511	//if(max(neighbours_wet[0], max(neighbours_wet[1],neighbours_wet[2]))>0){
512	// //usedouble == average edge stage
513	// usedouble = stage_edges[0]neighbours_wet[0] + stage_edges[1]neighbours_wet[1] + stage_edges[2]*neighbours_wet[2];
514	// usedouble /= 1.0*(neighbours_wet[0]+neighbours_wet[1]+neighbours_wet[2]);
515	//}else{
516	// usedouble = zl;
517	//}
518	//usedouble=min(usedouble,stage_centroid_values[k]);
519	//xmom_explicit_update[k] -= -0.5gmax(usedouble-zl,0.0)(usedouble-zl)normals[ki2]*length;
520	//ymom_explicit_update[k] -= -0.5gmax(usedouble-zl,0.0)(usedouble-zl)normals[ki2+1]*length;
521	////xmom_explicit_update[k] -= -0.5gabs(usedouble-zl)(usedouble-zl)normals[ki2]*length;
522	////ymom_explicit_update[k] -= -0.5gabs(usedouble-zl)(usedouble-zl)normals[ki2+1]*length;
523
524	//}
525
526	already_computed_flux[ki] = call; // #k Done
527
528
529	// Update neighbour n with same flux but reversed sign
530	if (n >= 0) {
531	stage_explicit_update[n] += edgeflux[0];
532	xmom_explicit_update[n] += edgeflux[1];
533	ymom_explicit_update[n] += edgeflux[2];
534	//Add bed slope term here
535	if((stage_centroid_values[n]>bed_centroid_values[n]+0.0)){
536	bedslope_work = -0.5gmax(stage_centroid_values[n]-zr,0.)(stage_centroid_values[n]-zr)length;
537	xmom_explicit_update[n] += normals[ki2]*bedslope_work;
538	ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
539	}else{
540	// Treat nearly dry cells
541	//bedslope_work = -0.5gmax(qr[0]-zr,0.)(qr[0]-zr)length; //-pressure_fluxlength; //-0.5gmax(qr[0]-zr,0.)(qr[0]-zr)*length;
542	bedslope_work = -pressure_flux*length;
543	xmom_explicit_update[n] += normals[ki2]*bedslope_work;
544	ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
545	}
546	// The triangle is (partially?) dry
547	// Find the average stage of the wet neighbours edges -- set this as the within triangle stage for the bedslope term
548	// FIXME: Note that we keep recomputing this -- could improve efficiency
549	//for(j=0;j<3;j++){
550	// // For each edge, check if each neighbouring triangle is wet
551	// useint = neighbours[3*n+j];
552	// if(useint>0){
553	// // neighbour is not a boundary condition
554	// if(stage_centroid_values[useint] > bed_centroid_values[useint]){
555	// //wet
556	// neighbours_wet[j]=1;
557	// useint = useint3 + neighbour_edges[3n+j]; // Index of neighbouring edge
558	// stage_edges[j] = stage_edge_values[useint]; // stage of neighbouring edge
559	// }else{
560	// //dry
561	// neighbours_wet[j]=0;
562	// stage_edges[j]=0.0; //Arbitrary
563	// }
564	// }else{
565	// // neighbour is a boundary condition, ASSUME IT IS WET
566	// useint = -useint -1;
567	// neighbours_wet[j]=1;
568	// stage_edges[j] = stage_boundary_values[useint];
569	// }
570	//}
571	////calculate average neighbouring edge stage, only considering wet triangles
572	//if(max(neighbours_wet[0],max(neighbours_wet[1],neighbours_wet[2]))>0){
573	// //average edge stage
574	// usedouble = stage_edges[0]neighbours_wet[0] + stage_edges[1]neighbours_wet[1] + stage_edges[2]*neighbours_wet[2];
575	// usedouble /= 1.0*(neighbours_wet[0]+neighbours_wet[1]+neighbours_wet[2]);
576	//}else{
577	// usedouble = zr;
578	//}
579	//usedouble=min(usedouble,stage_centroid_values[n]);
580	//xmom_explicit_update[n] += -0.5gmax(usedouble-zr,0.0)(usedouble-zr)normals[ki2]*length;
581	//ymom_explicit_update[n] += -0.5gmax(usedouble-zr,0.0)(usedouble-zr)normals[ki2+1]*length;
582	////xmom_explicit_update[n] += -0.5gabs(usedouble-zr)(usedouble-zr)normals[ki2]*length;
583	////ymom_explicit_update[n] += -0.5gabs(usedouble-zr)(usedouble-zr)normals[ki2+1]*length;
584
585	//}
586
587	already_computed_flux[nm] = call; // #n Done
588	}
589
590	// Update timestep based on edge i and possibly neighbour n
591	if (tri_full_flag[k] == 1) {
592	if (max_speed > epsilon) {
593	// Apply CFL condition for triangles joining this edge (triangle k and triangle n)
594
595	// CFL for triangle k
596	timestep = min(timestep, radii[k] / max_speed);
597
598	if (n >= 0) {
599	// Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
600	timestep = min(timestep, radii[n] / max_speed);
601	}
602
603	// Ted Rigby's suggested less conservative version
604	//if(n>=0){
605	// timestep = min(timestep, (radii[k]+radii[n])/max_speed);
606	//}else{
607	// timestep = min(timestep, radii[k]/max_speed);
608	//}
609	}
610	}
611
612	} // End edge i (and neighbour n)
613
614
615	// Normalise triangle k by area and store for when all conserved
616	// quantities get updated
617	inv_area = 1.0 / areas[k];
618	stage_explicit_update[k] *= inv_area;
619	xmom_explicit_update[k] *= inv_area;
620	ymom_explicit_update[k] *= inv_area;
621
622
623	// Keep track of maximal speeds
624	max_speed_array[k] = max_speed;
625
626	} // End triangle k
627
628
629	// Try to adjust stage_explicit_update to account for the negative stages,
630	// using an approach assumes a very small discharge is enough to allow
631	// the depth to increase from negative to zero
632	//for(k=0;k<number_of_elements;k++){
633	//depthtemp = stage_centroid_values[k]-bed_centroid_values[k];
634	// depthtemp = stage_centroid_values[k]-max(stage_edge_values[3k+1], max(stage_edge_values[3k+2], stage_edge_values[3*k]));
635	// if(depthtemp< 0.0){
636	// depthtemp = -depthtemp; // stage defecit
637	// usedouble = 2.0; //1.0e+2; //Large constant
638	// if(stage_explicit_update[k]timestepusedouble > depthtemp){
639	// stage_explicit_update[k] += (depthtemp - depthtemp/usedouble)/timestep;
640	// }else{
641	// stage_explicit_update[k] *=usedouble ;
642	// }
643	// }
644	//}
645
646	return timestep;
647	}
648
649	// Protect against the water elevation falling below the triangle bed
650	int _protect(int N,
651	double minimum_allowed_height,
652	double maximum_allowed_speed,
653	double epsilon,
654	double* wc,
655	double* zc,
656	double* zv,
657	double* xmomc,
658	double* ymomc) {
659
660	int k;
661	double hc, wmin;
662	double u, v, reduced_speed;
663
664
665	// Protect against initesimal and negative heights
666	if (maximum_allowed_speed < epsilon) {
667	for (k=0; k<N; k++) {
668	hc = wc[k] - zc[k];
669
670	if (hc < minimum_allowed_height) {
671
672	// Set momentum to zero and ensure h is non negative
673	xmomc[k] = 0.0;
674	ymomc[k] = 0.0;
675	if (hc <= 0.0){
676	// Definine the minimum possible stage value as the minimum bed edge value on triangle k.
677	wmin = 0.5min((zv[3k]+zv[3k+1]),min((zv[3k+1]+zv[3k+2]),(zv[3k+2]+zv[3*k])));
678	wc[k] = max(wc[k], wmin);
679	}
680	}
681	}
682	} else {
683
684	// Protect against initesimal and negative heights
685	for (k=0; k<N; k++) {
686	hc = wc[k] - zc[k];
687
688	if (hc < minimum_allowed_height) {
689
690	//New code: Adjust momentum to guarantee speeds are physical
691	// ensure h is non negative
692
693	if (hc <= 0.0) {
694	wc[k] = zc[k];
695	xmomc[k] = 0.0;
696	ymomc[k] = 0.0;
697	} else {
698	//Reduce excessive speeds derived from division by small hc
699	//FIXME (Ole): This may be unnecessary with new slope limiters
700	//in effect.
701
702	u = xmomc[k]/hc;
703	if (fabs(u) > maximum_allowed_speed) {
704	reduced_speed = maximum_allowed_speed * u/fabs(u);
705	//printf("Speed (u) has been reduced from %.3f to %.3f\n",
706	// u, reduced_speed);
707	xmomc[k] = reduced_speed * hc;
708	}
709
710	v = ymomc[k]/hc;
711	if (fabs(v) > maximum_allowed_speed) {
712	reduced_speed = maximum_allowed_speed * v/fabs(v);
713	//printf("Speed (v) has been reduced from %.3f to %.3f\n",
714	// v, reduced_speed);
715	ymomc[k] = reduced_speed * hc;
716	}
717	}
718	}
719	}
720	}
721	return 0;
722	}
723
724	//=========================================================================
725	// Python Glue
726	//=========================================================================
727
728
729	//========================================================================
730	// Compute fluxes
731	//========================================================================
732
733	// Modified central flux function
734
735	PyObject compute_fluxes_ext_central(PyObject self, PyObject *args) {
736	/*Compute all fluxes and the timestep suitable for all volumes
737	in domain.
738
739	Compute total flux for each conserved quantity using "flux_function_central"
740
741	Fluxes across each edge are scaled by edgelengths and summed up
742	Resulting flux is then scaled by area and stored in
743	explicit_update for each of the three conserved quantities
744	stage, xmomentum and ymomentum
745
746	The maximal allowable speed computed by the flux_function for each volume
747	is converted to a timestep that must not be exceeded. The minimum of
748	those is computed as the next overall timestep.
749
750	Python call:
751	domain.timestep = compute_fluxes(timestep,
752	domain.epsilon,
753	domain.H0,
754	domain.g,
755	domain.neighbours,
756	domain.neighbour_edges,
757	domain.normals,
758	domain.edgelengths,
759	domain.radii,
760	domain.areas,
761	tri_full_flag,
762	Stage.edge_values,
763	Xmom.edge_values,
764	Ymom.edge_values,
765	Bed.edge_values,
766	Stage.boundary_values,
767	Xmom.boundary_values,
768	Ymom.boundary_values,
769	Stage.explicit_update,
770	Xmom.explicit_update,
771	Ymom.explicit_update,
772	already_computed_flux,
773	optimise_dry_cells,
774	stage.centroid_values,
775	bed.centroid_values)
776
777
778	Post conditions:
779	domain.explicit_update is reset to computed flux values
780	domain.timestep is set to the largest step satisfying all volumes.
781
782
783	*/
784
785
786	PyArrayObject neighbours, neighbour_edges,
787	normals, edgelengths, radii, areas,
788	*tri_full_flag,
789	*stage_edge_values,
790	*xmom_edge_values,
791	*ymom_edge_values,
792	*bed_edge_values,
793	*stage_boundary_values,
794	*xmom_boundary_values,
795	*ymom_boundary_values,
796	*stage_explicit_update,
797	*xmom_explicit_update,
798	*ymom_explicit_update,
799	*already_computed_flux, //Tracks whether the flux across an edge has already been computed
800	*max_speed_array, //Keeps track of max speeds for each triangle
801	*stage_centroid_values,
802	*bed_centroid_values,
803	*bed_vertex_values;
804
805	double timestep, epsilon, H0, g;
806	int optimise_dry_cells;
807
808	// Convert Python arguments to C
809	if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOiOOO",
810	&timestep,
811	&epsilon,
812	&H0,
813	&g,
814	&neighbours,
815	&neighbour_edges,
816	&normals,
817	&edgelengths, &radii, &areas,
818	&tri_full_flag,
819	&stage_edge_values,
820	&xmom_edge_values,
821	&ymom_edge_values,
822	&bed_edge_values,
823	&stage_boundary_values,
824	&xmom_boundary_values,
825	&ymom_boundary_values,
826	&stage_explicit_update,
827	&xmom_explicit_update,
828	&ymom_explicit_update,
829	&already_computed_flux,
830	&max_speed_array,
831	&optimise_dry_cells,
832	&stage_centroid_values,
833	&bed_centroid_values,
834	&bed_vertex_values)) {
835	report_python_error(AT, "could not parse input arguments");
836	return NULL;
837	}
838
839	// check that numpy array objects arrays are C contiguous memory
840	CHECK_C_CONTIG(neighbours);
841	CHECK_C_CONTIG(neighbour_edges);
842	CHECK_C_CONTIG(normals);
843	CHECK_C_CONTIG(edgelengths);
844	CHECK_C_CONTIG(radii);
845	CHECK_C_CONTIG(areas);
846	CHECK_C_CONTIG(tri_full_flag);
847	CHECK_C_CONTIG(stage_edge_values);
848	CHECK_C_CONTIG(xmom_edge_values);
849	CHECK_C_CONTIG(ymom_edge_values);
850	CHECK_C_CONTIG(bed_edge_values);
851	CHECK_C_CONTIG(stage_boundary_values);
852	CHECK_C_CONTIG(xmom_boundary_values);
853	CHECK_C_CONTIG(ymom_boundary_values);
854	CHECK_C_CONTIG(stage_explicit_update);
855	CHECK_C_CONTIG(xmom_explicit_update);
856	CHECK_C_CONTIG(ymom_explicit_update);
857	CHECK_C_CONTIG(already_computed_flux);
858	CHECK_C_CONTIG(max_speed_array);
859	CHECK_C_CONTIG(stage_centroid_values);
860	CHECK_C_CONTIG(bed_centroid_values);
861	CHECK_C_CONTIG(bed_vertex_values);
862
863	int number_of_elements = stage_edge_values -> dimensions[0];
864
865	// Call underlying flux computation routine and update
866	// the explicit update arrays
867	timestep = _compute_fluxes_central(number_of_elements,
868	timestep,
869	epsilon,
870	H0,
871	g,
872	(long*) neighbours -> data,
873	(long*) neighbour_edges -> data,
874	(double*) normals -> data,
875	(double*) edgelengths -> data,
876	(double*) radii -> data,
877	(double*) areas -> data,
878	(long*) tri_full_flag -> data,
879	(double*) stage_edge_values -> data,
880	(double*) xmom_edge_values -> data,
881	(double*) ymom_edge_values -> data,
882	(double*) bed_edge_values -> data,
883	(double*) stage_boundary_values -> data,
884	(double*) xmom_boundary_values -> data,
885	(double*) ymom_boundary_values -> data,
886	(double*) stage_explicit_update -> data,
887	(double*) xmom_explicit_update -> data,
888	(double*) ymom_explicit_update -> data,
889	(long*) already_computed_flux -> data,
890	(double*) max_speed_array -> data,
891	optimise_dry_cells,
892	(double*) stage_centroid_values -> data,
893	(double*) bed_centroid_values -> data,
894	(double*) bed_vertex_values -> data);
895
896	// Return updated flux timestep
897	return Py_BuildValue("d", timestep);
898	}
899
900
901	PyObject flux_function_central(PyObject self, PyObject *args) {
902	//
903	// Gateway to innermost flux function.
904	// This is only used by the unit tests as the c implementation is
905	// normally called by compute_fluxes in this module.
906
907
908	PyArrayObject normal, ql, qr, edgeflux;
909	double g, epsilon, max_speed, H0, zl, zr;
910	double h0, limiting_threshold, pressure_flux;
911
912	if (!PyArg_ParseTuple(args, "OOOddOddd",
913	&normal, &ql, &qr, &zl, &zr, &edgeflux,
914	&epsilon, &g, &H0)) {
915	report_python_error(AT, "could not parse input arguments");
916	return NULL;
917	}
918
919
920	h0 = H0*H0; // This ensures a good balance when h approaches H0.
921	// But evidence suggests that h0 can be as little as
922	// epsilon!
923
924	limiting_threshold = 10*H0; // Avoid applying limiter below this
925	// threshold for performance reasons.
926	// See ANUGA manual under flux limiting
927
928	pressure_flux = 0.0; // Store the water pressure related component of the flux
929	_flux_function_central((double*) ql -> data,
930	(double*) qr -> data,
931	zl,
932	zr,
933	((double*) normal -> data)[0],
934	((double*) normal -> data)[1],
935	epsilon, h0, limiting_threshold,
936	g,
937	(double*) edgeflux -> data,
938	&max_speed,
939	&pressure_flux);
940
941	return Py_BuildValue("d", max_speed);
942	}
943
944	//========================================================================
945	// Gravity
946	//========================================================================
947
948	PyObject gravity(PyObject self, PyObject *args) {
949	//
950	// gravity(g, h, v, x, xmom, ymom)
951	//
952
953
954	PyArrayObject h, v, x, xmom, *ymom;
955	int k, N, k3, k6;
956	double g, avg_h, zx, zy;
957	double x0, y0, x1, y1, x2, y2, z0, z1, z2;
958	//double epsilon;
959
960	if (!PyArg_ParseTuple(args, "dOOOOO",
961	&g, &h, &v, &x,
962	&xmom, &ymom)) {
963	//&epsilon)) {
964	PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
965	return NULL;
966	}
967
968	// check that numpy array objects arrays are C contiguous memory
969	CHECK_C_CONTIG(h);
970	CHECK_C_CONTIG(v);
971	CHECK_C_CONTIG(x);
972	CHECK_C_CONTIG(xmom);
973	CHECK_C_CONTIG(ymom);
974
975	N = h -> dimensions[0];
976	for (k=0; k<N; k++) {
977	k3 = 3*k; // base index
978
979	// Get bathymetry
980	z0 = ((double*) v -> data)[k3 + 0];
981	z1 = ((double*) v -> data)[k3 + 1];
982	z2 = ((double*) v -> data)[k3 + 2];
983
984	// Optimise for flat bed
985	// Note (Ole): This didn't produce measurable speed up.
986	// Revisit later
987	//if (fabs(z0-z1)<epsilon && fabs(z1-z2)<epsilon) {
988	// continue;
989	//}
990
991	// Get average depth from centroid values
992	avg_h = ((double *) h -> data)[k];
993
994	// Compute bed slope
995	k6 = 6*k; // base index
996
997	x0 = ((double*) x -> data)[k6 + 0];
998	y0 = ((double*) x -> data)[k6 + 1];
999	x1 = ((double*) x -> data)[k6 + 2];
1000	y1 = ((double*) x -> data)[k6 + 3];
1001	x2 = ((double*) x -> data)[k6 + 4];
1002	y2 = ((double*) x -> data)[k6 + 5];
1003
1004
1005	_gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
1006
1007	// Update momentum
1008	((double) xmom -> data)[k] += -gzx*avg_h;
1009	((double) ymom -> data)[k] += -gzy*avg_h;
1010	}
1011
1012	return Py_BuildValue("");
1013	}
1014
1015
1016	//========================================================================
1017	// Protect -- to prevent the water level from falling below the minimum
1018	// bed_edge_value
1019	//========================================================================
1020
1021	PyObject protect(PyObject self, PyObject *args) {
1022	//
1023	// protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
1024
1025
1026	PyArrayObject
1027	*wc, //Stage at centroids
1028	*zc, //Elevation at centroids
1029	*zv, //Elevation at vertices
1030	*xmomc, //Momentums at centroids
1031	*ymomc;
1032
1033
1034	int N;
1035	double minimum_allowed_height, maximum_allowed_speed, epsilon;
1036
1037	// Convert Python arguments to C
1038	if (!PyArg_ParseTuple(args, "dddOOOOO",
1039	&minimum_allowed_height,
1040	&maximum_allowed_speed,
1041	&epsilon,
1042	&wc, &zc, &zv, &xmomc, &ymomc)) {
1043	report_python_error(AT, "could not parse input arguments");
1044	return NULL;
1045	}
1046
1047	N = wc -> dimensions[0];
1048
1049	_protect(N,
1050	minimum_allowed_height,
1051	maximum_allowed_speed,
1052	epsilon,
1053	(double*) wc -> data,
1054	(double*) zc -> data,
1055	(double*) zv -> data,
1056	(double*) xmomc -> data,
1057	(double*) ymomc -> data);
1058
1059	return Py_BuildValue("");
1060	}
1061
1062	//========================================================================
1063	// Method table for python module
1064	//========================================================================
1065
1066	static struct PyMethodDef MethodTable[] = {
1067	/* The cast of the function is necessary since PyCFunction values
1068	* only take two PyObject* parameters, and rotate() takes
1069	* three.
1070	*/
1071	{"compute_fluxes_ext_central", compute_fluxes_ext_central, METH_VARARGS, "Print out"},
1072	{"gravity_c", gravity, METH_VARARGS, "Print out"},
1073	{"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},
1074	{"protect", protect, METH_VARARGS \| METH_KEYWORDS, "Print out"},
1075	{NULL, NULL, 0, NULL}
1076	};
1077
1078	// Module initialisation
1079	void initswb2_domain_ext(void){
1080	Py_InitModule("swb2_domain_ext", MethodTable);
1081
1082	import_array(); // Necessary for handling of NumPY structures
1083	}

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

Download in other formats: