Context Navigation

source: anuga_work/development/anuga_1d/comp_flux_vel_ext.c @ 7815

Last change on this file since 7815 was 6042, checked in by steve, 16 years ago
Found bug associated with changing max_speed_array to an integer
File size: 17.7 KB

Line
1	#include "Python.h"
2	#include "Numeric/arrayobject.h"
3	#include "math.h"
4	#include <stdio.h>
5	const double pi = 3.14159265358979;
6
7
8	// Shared code snippets
9	#include "util_ext.h"
10
11
12	/* double max(double a, double b) { */
13	/* double z; */
14	/* z=(a>b)?a:b; */
15	/* return z;} */
16
17	/* double min(double a, double b) { */
18	/* double z; */
19	/* z=(a<b)?a:b; */
20	/* return z;} */
21
22
23	// Function to obtain speed from momentum and depth.
24	// This is used by flux functions
25	// Input parameters uh and h may be modified by this function.
26	double _compute_speed(double *uh,
27	double *h,
28	double epsilon,
29	double h0) {
30
31	double u;
32
33	if (*h < epsilon) {
34	*h = 0.0; //Could have been negative
35	u = 0.0;
36	} else {
37	u = uh/(h + h0/ *h);
38	}
39
40
41	// Adjust momentum to be consistent with speed
42	uh = u *h;
43
44	return u;
45	}
46
47	//-------------------------------------------------------------
48	// New vel based code
49	//-------------------------------------------------------------
50	//Innermost flux function (using w=z+h)
51	int _flux_function_vel(double q_left, double q_right,
52	double normals, double g, double epsilon, double h0,
53	double edgeflux, double max_speed) {
54
55	int i;
56	double flux_left[2], flux_right[2];
57	double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;
58	double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;
59	double z, s_max, s_min, denom;
60
61
62	w_left = q_left[0];
63	uh_left = q_left[1]*normals;
64	z_left = q_left[2];
65	h_left = q_left[3];
66	u_left = q_left[4]*normals;
67
68	/* printf("w_left = %f \n",w_left); */
69	/* printf("uh_left = %f \n",uh_left); */
70	/* printf("z_left = %f \n",z_left); */
71	/* printf("h_left = %f \n",h_left); */
72	/* printf("u_left = %f \n",u_left); */
73
74	w_right = q_right[0];
75	uh_right = q_right[1]*normals;
76	z_right = q_right[2];
77	h_right = q_right[3];
78	u_right = q_right[4]*normals;
79
80	z = (z_left+z_right)/2.0;
81
82	soundspeed_left = sqrt(g*h_left);
83	soundspeed_right = sqrt(g*h_right);
84
85	s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
86	if (s_max < 0.0) s_max = 0.0;
87
88	s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
89	if (s_min > 0.0) s_min = 0.0;
90
91
92	// Flux formulas
93	flux_left[0] = u_left*h_left;
94	flux_left[1] = u_leftu_lefth_left + 0.5gh_left*h_left;
95
96	flux_right[0] = u_right*h_right;
97	flux_right[1] = u_rightu_righth_right + 0.5gh_right*h_right;
98
99	// Flux computation
100	denom = s_max-s_min;
101	if (denom < epsilon) {
102	for (i=0; i<2; i++) edgeflux[i] = 0.0;
103	*max_speed = 0.0;
104	} else {
105	edgeflux[0] = s_maxflux_left[0] - s_minflux_right[0];
106	edgeflux[0] += s_maxs_min(w_right-w_left);
107	edgeflux[0] /= denom;
108	edgeflux[1] = s_maxflux_left[1] - s_minflux_right[1];
109	edgeflux[1] += s_maxs_min(uh_right-uh_left);
110	edgeflux[1] /= denom;
111	edgeflux[1] *= normals;
112
113	// Maximal wavespeed
114	*max_speed = max(fabs(s_max), fabs(s_min));
115	}
116	return 0;
117	}
118
119	// Computational function for flux computation
120	double _compute_fluxes_vel_ext(double cfl,
121	double timestep,
122	double epsilon,
123	double g,
124	double h0,
125	long* neighbours,
126	long* neighbour_vertices,
127	double* normals,
128	double* areas,
129	double* stage_edge_values,
130	double* xmom_edge_values,
131	double* bed_edge_values,
132	double* height_edge_values,
133	double* velocity_edge_values,
134	double* stage_boundary_values,
135	double* xmom_boundary_values,
136	double* bed_boundary_values,
137	double* height_boundary_values,
138	double* velocity_boundary_values,
139	double* stage_explicit_update,
140	double* xmom_explicit_update,
141	int number_of_elements,
142	double* max_speed_array) {
143
144	double flux[2], ql[5], qr[5], edgeflux[2];
145	double max_speed, normal;
146	int k, i, ki, n, m, nm=0;
147
148	//printf("Inside _comp\n");
149
150	for (k=0; k<number_of_elements; k++) {
151	flux[0] = 0.0;
152	flux[1] = 0.0;
153	//printf("k = %d\n",k);
154
155	for (i=0; i<2; i++) {
156	ki = k*2+i;
157
158	ql[0] = stage_edge_values[ki];
159	ql[1] = xmom_edge_values[ki];
160	ql[2] = bed_edge_values[ki];
161	ql[3] = height_edge_values[ki];
162	ql[4] = velocity_edge_values[ki];
163
164	n = neighbours[ki];
165	if (n<0) {
166	m = -n-1;
167	qr[0] = stage_boundary_values[m];
168	qr[1] = xmom_boundary_values[m];
169	qr[2] = bed_boundary_values[m];
170	qr[3] = height_boundary_values[m];
171	qr[4] = velocity_boundary_values[m];
172	} else {
173	m = neighbour_vertices[ki];
174	nm = n*2+m;
175	qr[0] = stage_edge_values[nm];
176	qr[1] = xmom_edge_values[nm];
177	qr[2] = bed_edge_values[nm];
178	qr[3] = height_edge_values[nm];
179	qr[4] = velocity_edge_values[nm];
180	}
181
182	normal = normals[ki];
183	_flux_function_vel(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);
184	flux[0] -= edgeflux[0];
185	flux[1] -= edgeflux[1];
186
187	// Update timestep based on edge i and possibly neighbour n
188	if (max_speed > epsilon) {
189	// Original CFL calculation
190
191	timestep = min(timestep, 0.5cflareas[k]/max_speed);
192	if (n>=0) {
193	timestep = min(timestep, 0.5cflareas[n]/max_speed);
194	}
195	}
196	} // End edge i (and neighbour n)
197	flux[0] /= areas[k];
198	stage_explicit_update[k] = flux[0];
199	flux[1] /= areas[k];
200	xmom_explicit_update[k] = flux[1];
201
202	//Keep track of maximal speeds
203	max_speed_array[k]=max_speed;
204	}
205	return timestep;
206	}
207
208	//-------------------------------------------------------------
209	// Old code
210	//------------------------------------------------------------
211	//Innermost flux function (using w=z+h)
212	int _flux_function(double q_left, double q_right,
213	double z_left, double z_right,
214	double normals, double g, double epsilon, double h0,
215	double edgeflux, double max_speed) {
216
217	int i;
218	double ql[2], qr[2], flux_left[2], flux_right[2];
219	double z, w_left, h_left, uh_left, soundspeed_left, u_left;
220	double w_right, h_right, uh_right, soundspeed_right, u_right;
221	double s_max, s_min, denom;
222
223	//printf("h0 = %f \n",h0);
224	ql[0] = q_left[0];
225	ql[1] = q_left[1];
226	ql[1] = ql[1]*normals;
227
228	qr[0] = q_right[0];
229	qr[1] = q_right[1];
230	qr[1] = qr[1]*normals;
231
232	z = (z_left+z_right)/2.0;
233
234	//w_left = ql[0];
235	//h_left = w_left-z;
236	//uh_left = ql[1];
237
238
239
240	// Compute speeds in x-direction
241	w_left = ql[0];
242	h_left = w_left-z;
243	uh_left = ql[1];
244
245	u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
246
247	w_right = qr[0];
248	h_right = w_right-z;
249	uh_right = qr[1];
250
251	u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
252
253	soundspeed_left = sqrt(g*h_left);
254	soundspeed_right = sqrt(g*h_right);
255
256	s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
257	if (s_max < 0.0) s_max = 0.0;
258
259	s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
260	if (s_min > 0.0) s_min = 0.0;
261
262
263	// Flux formulas
264	flux_left[0] = u_left*h_left;
265	flux_left[1] = u_leftuh_left + 0.5gh_lefth_left;
266
267	flux_right[0] = u_right*h_right;
268	flux_right[1] = u_rightuh_right + 0.5gh_righth_right;
269
270	// Flux computation
271	denom = s_max-s_min;
272	if (denom < epsilon) {
273	for (i=0; i<2; i++) edgeflux[i] = 0.0;
274	*max_speed = 0.0;
275	} else {
276	edgeflux[0] = s_maxflux_left[0] - s_minflux_right[0];
277	edgeflux[0] += s_maxs_min(qr[0]-ql[0]);
278	edgeflux[0] /= denom;
279	edgeflux[1] = s_maxflux_left[1] - s_minflux_right[1];
280	edgeflux[1] += s_maxs_min(qr[1]-ql[1]);
281	edgeflux[1] /= denom;
282	edgeflux[1] *= normals;
283
284	// Maximal wavespeed
285	*max_speed = max(fabs(s_max), fabs(s_min));
286	}
287	return 0;
288	}
289
290
291
292
293	// Computational function for flux computation
294	double _compute_fluxes_ext(
295	double cfl,
296	double timestep,
297	double epsilon,
298	double g,
299	double h0,
300	long* neighbours,
301	long* neighbour_vertices,
302	double* normals,
303	double* areas,
304	double* stage_edge_values,
305	double* xmom_edge_values,
306	double* bed_edge_values,
307	double* stage_boundary_values,
308	double* xmom_boundary_values,
309	double* stage_explicit_update,
310	double* xmom_explicit_update,
311	int number_of_elements,
312	double* max_speed_array) {
313
314	double flux[2], ql[2], qr[2], edgeflux[2];
315	double zl, zr, max_speed, normal;
316	int k, i, ki, n, m, nm=0;
317
318
319	for (k=0; k<number_of_elements; k++) {
320	flux[0] = 0.0;
321	flux[1] = 0.0;
322
323	for (i=0; i<2; i++) {
324	ki = k*2+i;
325
326	ql[0] = stage_edge_values[ki];
327	ql[1] = xmom_edge_values[ki];
328	zl = bed_edge_values[ki];
329
330	n = neighbours[ki];
331	if (n<0) {
332	m = -n-1;
333	qr[0] = stage_boundary_values[m];
334	qr[1] = xmom_boundary_values[m];
335	zr = zl;
336	} else {
337	m = neighbour_vertices[ki];
338	nm = n*2+m;
339	qr[0] = stage_edge_values[nm];
340	qr[1] = xmom_edge_values[nm];
341	zr = bed_edge_values[nm];
342	}
343
344	normal = normals[ki];
345	_flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
346	flux[0] -= edgeflux[0];
347	flux[1] -= edgeflux[1];
348
349	// Update timestep based on edge i and possibly neighbour n
350	if (max_speed > epsilon) {
351	// Original CFL calculation
352
353	timestep = min(timestep, 0.5cflareas[k]/max_speed);
354	if (n>=0) {
355	timestep = min(timestep, 0.5cflareas[n]/max_speed);
356	}
357	}
358	} // End edge i (and neighbour n)
359	flux[0] /= areas[k];
360	stage_explicit_update[k] = flux[0];
361	flux[1] /= areas[k];
362	xmom_explicit_update[k] = flux[1];
363
364	//Keep track of maximal speeds
365	max_speed_array[k]=max_speed;
366	}
367	return timestep;
368	}
369
370	//=========================================================================
371	// Python Glue
372	//=========================================================================
373	PyObject compute_fluxes_ext(PyObject self, PyObject *args) {
374
375	PyObject
376	*domain,
377	*stage,
378	*xmom,
379	*bed;
380
381	PyArrayObject
382	*neighbours,
383	*neighbour_vertices,
384	*normals,
385	*areas,
386	*stage_vertex_values,
387	*xmom_vertex_values,
388	*bed_vertex_values,
389	*stage_boundary_values,
390	*xmom_boundary_values,
391	*stage_explicit_update,
392	*xmom_explicit_update,
393	*max_speed_array;
394
395	double timestep, epsilon, g, h0, cfl;
396	int number_of_elements;
397
398
399	// Convert Python arguments to C
400	if (!PyArg_ParseTuple(args, "dOOOO",
401	&timestep,
402	&domain,
403	&stage,
404	&xmom,
405	&bed)) {
406	PyErr_SetString(PyExc_RuntimeError, "comp_flux_vel_ext.c: compute_fluxes_ext could not parse input");
407	return NULL;
408	}
409
410
411	epsilon = get_python_double(domain,"epsilon");
412	g = get_python_double(domain,"g");
413	h0 = get_python_double(domain,"h0");
414	cfl = get_python_double(domain,"CFL");
415
416
417	neighbours = get_consecutive_array(domain, "neighbours");
418	neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
419	normals = get_consecutive_array(domain, "normals");
420	areas = get_consecutive_array(domain, "areas");
421	max_speed_array = get_consecutive_array(domain, "max_speed_array");
422
423	stage_vertex_values = get_consecutive_array(stage, "vertex_values");
424	xmom_vertex_values = get_consecutive_array(xmom, "vertex_values");
425	bed_vertex_values = get_consecutive_array(bed, "vertex_values");
426
427	stage_boundary_values = get_consecutive_array(stage, "boundary_values");
428	xmom_boundary_values = get_consecutive_array(xmom, "boundary_values");
429
430
431	stage_explicit_update = get_consecutive_array(stage, "explicit_update");
432	xmom_explicit_update = get_consecutive_array(xmom, "explicit_update");
433
434
435
436	number_of_elements = stage_vertex_values -> dimensions[0];
437
438
439
440	// Call underlying flux computation routine and update
441	// the explicit update arrays
442	timestep = _compute_fluxes_ext(
443	cfl,
444	timestep,
445	epsilon,
446	g,
447	h0,
448	(long*) neighbours -> data,
449	(long*) neighbour_vertices -> data,
450	(double*) normals -> data,
451	(double*) areas -> data,
452	(double*) stage_vertex_values -> data,
453	(double*) xmom_vertex_values -> data,
454	(double*) bed_vertex_values -> data,
455	(double*) stage_boundary_values -> data,
456	(double*) xmom_boundary_values -> data,
457	(double*) stage_explicit_update -> data,
458	(double*) xmom_explicit_update -> data,
459	number_of_elements,
460	(double*) max_speed_array -> data);
461
462
463	Py_DECREF(neighbours);
464	Py_DECREF(neighbour_vertices);
465	Py_DECREF(normals);
466	Py_DECREF(areas);
467	Py_DECREF(stage_vertex_values);
468	Py_DECREF(xmom_vertex_values);
469	Py_DECREF(bed_vertex_values);
470	Py_DECREF(stage_boundary_values);
471	Py_DECREF(xmom_boundary_values);
472	Py_DECREF(stage_explicit_update);
473	Py_DECREF(xmom_explicit_update);
474	Py_DECREF(max_speed_array);
475
476
477
478
479	// Return updated flux timestep
480	return Py_BuildValue("d", timestep);
481	}
482
483
484	//------------------------------------------------
485	// New velocity based compute fluxes
486	//------------------------------------------------
487	PyObject compute_fluxes_vel_ext(PyObject self, PyObject *args) {
488
489	PyObject
490	*domain,
491	*stage,
492	*xmom,
493	*bed,
494	*height,
495	*velocity;
496
497	PyArrayObject
498	*neighbours,
499	*neighbour_vertices,
500	*normals,
501	*areas,
502	*stage_vertex_values,
503	*xmom_vertex_values,
504	*bed_vertex_values,
505	*height_vertex_values,
506	*velocity_vertex_values,
507	*stage_boundary_values,
508	*xmom_boundary_values,
509	*bed_boundary_values,
510	*height_boundary_values,
511	*velocity_boundary_values,
512	*stage_explicit_update,
513	*xmom_explicit_update,
514	*max_speed_array;
515
516	double timestep, epsilon, g, h0, cfl;
517	int number_of_elements;
518
519
520	// Convert Python arguments to C
521	if (!PyArg_ParseTuple(args, "dOOOOOO",
522	&timestep,
523	&domain,
524	&stage,
525	&xmom,
526	&bed,
527	&height,
528	&velocity)) {
529	PyErr_SetString(PyExc_RuntimeError, "comp_flux_vel_ext.c: compute_fluxes_vel_ext could not parse input");
530	printf("comp_flux_vel_ext.c: compute_fluxes_vel_ext could not parse input");
531	return NULL;
532	}
533
534
535	epsilon = get_python_double(domain,"epsilon");
536	g = get_python_double(domain,"g");
537	h0 = get_python_double(domain,"h0");
538	cfl = get_python_double(domain,"CFL");
539
540
541	neighbours = get_consecutive_array(domain, "neighbours");
542	neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
543	normals = get_consecutive_array(domain, "normals");
544	areas = get_consecutive_array(domain, "areas");
545	max_speed_array = get_consecutive_array(domain, "max_speed_array");
546
547	stage_vertex_values = get_consecutive_array(stage, "vertex_values");
548	xmom_vertex_values = get_consecutive_array(xmom, "vertex_values");
549	bed_vertex_values = get_consecutive_array(bed, "vertex_values");
550	height_vertex_values = get_consecutive_array(height, "vertex_values");
551	velocity_vertex_values = get_consecutive_array(velocity, "vertex_values");
552
553	stage_boundary_values = get_consecutive_array(stage, "boundary_values");
554	xmom_boundary_values = get_consecutive_array(xmom, "boundary_values");
555	bed_boundary_values = get_consecutive_array(bed, "boundary_values");
556	height_boundary_values = get_consecutive_array(height, "boundary_values");
557	velocity_boundary_values = get_consecutive_array(velocity, "boundary_values");
558
559
560	stage_explicit_update = get_consecutive_array(stage, "explicit_update");
561	xmom_explicit_update = get_consecutive_array(xmom, "explicit_update");
562
563	number_of_elements = stage_vertex_values -> dimensions[0];
564
565
566	//printf("comp_flux_vel 1 N = %d %15.5e %15.5e %15.5e %15.5e %15.5e \n",number_of_elements,cfl,timestep,epsilon,g,h0);
567
568
569
570	// Call underlying flux computation routine and update
571	// the explicit update arrays
572	timestep = _compute_fluxes_vel_ext(cfl,
573	timestep,
574	epsilon,
575	g,
576	h0,
577	(long*) neighbours -> data,
578	(long*) neighbour_vertices -> data,
579	(double*) normals -> data,
580	(double*) areas -> data,
581	(double*) stage_vertex_values -> data,
582	(double*) xmom_vertex_values -> data,
583	(double*) bed_vertex_values -> data,
584	(double*) height_vertex_values -> data,
585	(double*) velocity_vertex_values -> data,
586	(double*) stage_boundary_values -> data,
587	(double*) xmom_boundary_values -> data,
588	(double*) bed_boundary_values -> data,
589	(double*) height_boundary_values -> data,
590	(double*) velocity_boundary_values -> data,
591	(double*) stage_explicit_update -> data,
592	(double*) xmom_explicit_update -> data,
593	number_of_elements,
594	(double*) max_speed_array -> data);
595
596	//printf("comp_flux_vel 2 \n");
597
598	Py_DECREF(neighbours);
599	Py_DECREF(neighbour_vertices);
600	Py_DECREF(normals);
601	Py_DECREF(areas);
602	Py_DECREF(stage_vertex_values);
603	Py_DECREF(xmom_vertex_values);
604	Py_DECREF(bed_vertex_values);
605	Py_DECREF(height_vertex_values);
606	Py_DECREF(velocity_vertex_values);
607	Py_DECREF(stage_boundary_values);
608	Py_DECREF(xmom_boundary_values);
609	Py_DECREF(bed_boundary_values);
610	Py_DECREF(height_boundary_values);
611	Py_DECREF(velocity_boundary_values);
612	Py_DECREF(stage_explicit_update);
613	Py_DECREF(xmom_explicit_update);
614	Py_DECREF(max_speed_array);
615
616	//printf("comp_flux_vel 3 \n");
617
618	// Return updated flux timestep
619	return Py_BuildValue("d", timestep);
620	}
621
622
623
624	//-------------------------------
625	// Method table for python module
626	//-------------------------------
627
628	static struct PyMethodDef MethodTable[] = {
629	{"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
630	{"compute_fluxes_vel_ext", compute_fluxes_vel_ext, METH_VARARGS, "Print out"},
631	{NULL, NULL}
632	};
633
634	// Module initialisation
635	void initcomp_flux_vel_ext(void){
636	Py_InitModule("comp_flux_vel_ext", MethodTable);
637	import_array();
638	}

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

Download in other formats: