Context Navigation

source: anuga_work/development/2010-projects/anuga_1d/sww/sww_vel_comp_flux_ext.c @ 7840

Last change on this file since 7840 was 7839, checked in by steve, 14 years ago
Changing name of 1d projects so that it will be easy to moveto the trunk
File size: 17.7 KB

Line
1	#include "Python.h"
2	#include "numpy/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 initsww_vel_comp_flux_ext(void){
636	Py_InitModule("sww_vel_comp_flux_ext", MethodTable);
637	import_array();
638	}

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

Download in other formats: