Context Navigation

source: anuga_work/development/anuga_1d/comp_flux_ext.c @ 5724

Last change on this file since 5724 was 5724, checked in by sudi, 15 years ago

File size: 10.2 KB

Rev	Line
[5563]	1	#include "Python.h"
	2	#include "Numeric/arrayobject.h"
	3	#include "math.h"
	4	#include <stdio.h>
	5	const double pi = 3.14159265358979;
	6
[5564]	7
[5563]	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
[5724]	22	<<<<<<< .mine
[5587]	23
	24
[5724]	25
	26
	27	=======
	28
	29
[5587]	30	// Function to obtain speed from momentum and depth.
	31	// This is used by flux functions
	32	// Input parameters uh and h may be modified by this function.
	33	double _compute_speed(double *uh,
	34	double *h,
	35	double epsilon,
	36	double h0) {
	37
	38	double u;
	39
	40	if (*h < epsilon) {
	41	*h = 0.0; //Could have been negative
	42	u = 0.0;
	43	} else {
	44	u = uh/(h + h0/ *h);
	45	}
	46
	47
	48	// Adjust momentum to be consistent with speed
	49	uh = u *h;
	50
	51	return u;
	52	}
	53
	54
	55
[5724]	56	>>>>>>> .r5723
[5563]	57	//Innermost flux function (using w=z+h)
	58	int _flux_function(double q_left, double q_right,
	59	double z_left, double z_right,
[5587]	60	double normals, double g, double epsilon, double h0,
[5563]	61	double edgeflux, double max_speed) {
	62
	63	int i;
	64	double ql[2], qr[2], flux_left[2], flux_right[2];
	65	double z, w_left, h_left, uh_left, soundspeed_left, u_left;
	66	double w_right, h_right, uh_right, soundspeed_right, u_right;
	67	double s_max, s_min, denom;
	68
[5587]	69	//printf("h0 = %f \n",h0);
[5563]	70	ql[0] = q_left[0];
	71	ql[1] = q_left[1];
	72	ql[1] = ql[1]*normals;
	73
	74	qr[0] = q_right[0];
	75	qr[1] = q_right[1];
	76	qr[1] = qr[1]*normals;
	77
	78	z = (z_left+z_right)/2.0;
	79
[5724]	80	//w_left = ql[0];
	81	//h_left = w_left-z;
	82	//uh_left = ql[1];
[5563]	83
	84
	85
	86	// Compute speeds in x-direction
	87	w_left = ql[0];
	88	h_left = w_left-z;
	89	uh_left = ql[1];
[5587]	90
	91	u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
	92
[5563]	93	w_right = qr[0];
	94	h_right = w_right-z;
	95	uh_right = qr[1];
[5587]	96
	97	u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
[5563]	98
	99	soundspeed_left = sqrt(g*h_left);
	100	soundspeed_right = sqrt(g*h_right);
	101
	102	s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
	103	if (s_max < 0.0) s_max = 0.0;
	104
	105	s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
	106	if (s_min > 0.0) s_min = 0.0;
	107
	108
	109	// Flux formulas
	110	flux_left[0] = u_left*h_left;
	111	flux_left[1] = u_leftuh_left + 0.5gh_lefth_left;
	112
	113	flux_right[0] = u_right*h_right;
	114	flux_right[1] = u_rightuh_right + 0.5gh_righth_right;
	115
	116	// Flux computation
	117	denom = s_max-s_min;
	118	if (denom < epsilon) {
	119	for (i=0; i<2; i++) edgeflux[i] = 0.0;
	120	*max_speed = 0.0;
	121	} else {
	122	edgeflux[0] = s_maxflux_left[0] - s_minflux_right[0];
	123	edgeflux[0] += s_maxs_min(qr[0]-ql[0]);
	124	edgeflux[0] /= denom;
	125	edgeflux[1] = s_maxflux_left[1] - s_minflux_right[1];
	126	edgeflux[1] += s_maxs_min(qr[1]-ql[1]);
	127	edgeflux[1] /= denom;
	128	edgeflux[1] *= normals;
	129
	130	// Maximal wavespeed
	131	*max_speed = max(fabs(s_max), fabs(s_min));
	132	}
	133	return 0;
	134	}
	135
	136
	137
	138
	139	// Computational function for flux computation
	140	double _compute_fluxes_ext(double timestep,
[5587]	141	double epsilon,
	142	double g,
	143	double h0,
	144	long* neighbours,
	145	long* neighbour_vertices,
	146	double* normals,
	147	double* areas,
	148	double* stage_edge_values,
	149	double* xmom_edge_values,
	150	double* bed_edge_values,
	151	double* stage_boundary_values,
	152	double* xmom_boundary_values,
	153	double* stage_explicit_update,
	154	double* xmom_explicit_update,
	155	int number_of_elements,
	156	double* max_speed_array) {
[5563]	157
[5587]	158	double flux[2], ql[2], qr[2], edgeflux[2];
[5563]	159	double zl, zr, max_speed, normal;
	160	int k, i, ki, n, m, nm=0;
	161
	162	for (k=0; k<number_of_elements; k++) {
	163	flux[0] = 0.0;
	164	flux[1] = 0.0;
	165
	166	for (i=0; i<2; i++) {
	167	ki = k*2+i;
	168
	169	ql[0] = stage_edge_values[ki];
	170	ql[1] = xmom_edge_values[ki];
	171	zl = bed_edge_values[ki];
	172
	173	n = neighbours[ki];
	174	if (n<0) {
	175	m = -n-1;
	176	qr[0] = stage_boundary_values[m];
	177	qr[1] = xmom_boundary_values[m];
	178	zr = zl;
	179	} else {
	180	m = neighbour_vertices[ki];
	181	nm = n*2+m;
	182	qr[0] = stage_edge_values[nm];
	183	qr[1] = xmom_edge_values[nm];
	184	zr = bed_edge_values[nm];
	185	}
	186
	187	normal = normals[ki];
[5587]	188	_flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
[5563]	189	flux[0] -= edgeflux[0];
	190	flux[1] -= edgeflux[1];
	191
	192	// Update timestep based on edge i and possibly neighbour n
	193	if (max_speed > epsilon) {
	194	// Original CFL calculation
	195
	196	timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
	197	if (n>=0) {
	198	timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
	199	}
	200	}
	201	} // End edge i (and neighbour n)
	202	flux[0] /= areas[k];
	203	stage_explicit_update[k] = flux[0];
	204	flux[1] /= areas[k];
	205	xmom_explicit_update[k] = flux[1];
	206
	207	//Keep track of maximal speeds
	208	max_speed_array[k]=max_speed;
	209	}
[5724]	210	return timestep;
	211	}
[5563]	212
	213
	214
	215
	216
	217
	218
	219
	220	//=========================================================================
	221	// Python Glue
	222	//=========================================================================
	223	PyObject compute_fluxes_ext(PyObject self, PyObject *args) {
	224
	225	PyArrayObject
	226	*neighbours,
	227	*neighbour_vertices,
	228	*normals,
	229	*areas,
	230	*stage_edge_values,
	231	*xmom_edge_values,
	232	*bed_edge_values,
	233	*stage_boundary_values,
	234	*xmom_boundary_values,
	235	*stage_explicit_update,
	236	*xmom_explicit_update,
	237	*max_speed_array;
	238
[5587]	239	double timestep, epsilon, g, h0;
[5563]	240	int number_of_elements;
	241
	242	// Convert Python arguments to C
[5587]	243	if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOiO",
[5563]	244	&timestep,
	245	&epsilon,
	246	&g,
[5587]	247	&h0,
[5563]	248	&neighbours,
	249	&neighbour_vertices,
	250	&normals,
	251	&areas,
	252	&stage_edge_values,
	253	&xmom_edge_values,
	254	&bed_edge_values,
	255	&stage_boundary_values,
	256	&xmom_boundary_values,
	257	&stage_explicit_update,
	258	&xmom_explicit_update,
	259	&number_of_elements,
	260	&max_speed_array)) {
	261	PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
	262	return NULL;
	263	}
	264
	265
	266	// Call underlying flux computation routine and update
	267	// the explicit update arrays
	268	timestep = _compute_fluxes_ext(timestep,
	269	epsilon,
	270	g,
[5587]	271	h0,
[5563]	272	(long*) neighbours -> data,
	273	(long*) neighbour_vertices -> data,
	274	(double*) normals -> data,
	275	(double*) areas -> data,
	276	(double*) stage_edge_values -> data,
	277	(double*) xmom_edge_values -> data,
	278	(double*) bed_edge_values -> data,
	279	(double*) stage_boundary_values -> data,
	280	(double*) xmom_boundary_values -> data,
	281	(double*) stage_explicit_update -> data,
	282	(double*) xmom_explicit_update -> data,
	283	number_of_elements,
	284	(double*) max_speed_array -> data);
	285
	286	// Return updated flux timestep
	287	return Py_BuildValue("d", timestep);
	288	}
	289
	290
	291	PyObject compute_fluxes_ext_short(PyObject self, PyObject *args) {
	292
	293	PyObject
	294	*domain,
	295	*stage,
	296	*xmom,
	297	*bed;
	298
	299	PyArrayObject
	300	*neighbours,
	301	*neighbour_vertices,
	302	*normals,
	303	*areas,
	304	*stage_vertex_values,
	305	*xmom_vertex_values,
	306	*bed_vertex_values,
	307	*stage_boundary_values,
	308	*xmom_boundary_values,
	309	*stage_explicit_update,
	310	*xmom_explicit_update,
	311	*max_speed_array;
	312
[5587]	313	double timestep, epsilon, g, h0;
[5563]	314	int number_of_elements;
	315
	316
	317	// Convert Python arguments to C
	318	if (!PyArg_ParseTuple(args, "dOOOO",
	319	&timestep,
	320	&domain,
	321	&stage,
	322	&xmom,
	323	&bed)) {
	324	PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
	325	return NULL;
	326	}
	327
	328
	329	epsilon = get_python_double(domain,"epsilon");
	330	g = get_python_double(domain,"g");
[5587]	331	h0 = get_python_double(domain,"h0");
[5563]	332
	333
	334	neighbours = get_consecutive_array(domain, "neighbours");
	335	neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
	336	normals = get_consecutive_array(domain, "normals");
	337	areas = get_consecutive_array(domain, "areas");
	338	max_speed_array = get_consecutive_array(domain, "max_speed_array");
	339
	340	stage_vertex_values = get_consecutive_array(stage, "vertex_values");
	341	xmom_vertex_values = get_consecutive_array(xmom, "vertex_values");
	342	bed_vertex_values = get_consecutive_array(bed, "vertex_values");
	343
	344	stage_boundary_values = get_consecutive_array(stage, "boundary_values");
	345	xmom_boundary_values = get_consecutive_array(xmom, "boundary_values");
	346
	347
	348	stage_explicit_update = get_consecutive_array(stage, "explicit_update");
	349	xmom_explicit_update = get_consecutive_array(xmom, "explicit_update");
	350
	351
	352
	353	number_of_elements = stage_vertex_values -> dimensions[0];
	354
	355
	356
	357	// Call underlying flux computation routine and update
	358	// the explicit update arrays
	359	timestep = _compute_fluxes_ext(timestep,
[5587]	360	epsilon,
	361	g,
	362	h0,
	363	(long*) neighbours -> data,
	364	(long*) neighbour_vertices -> data,
	365	(double*) normals -> data,
	366	(double*) areas -> data,
	367	(double*) stage_vertex_values -> data,
	368	(double*) xmom_vertex_values -> data,
	369	(double*) bed_vertex_values -> data,
	370	(double*) stage_boundary_values -> data,
	371	(double*) xmom_boundary_values -> data,
	372	(double*) stage_explicit_update -> data,
	373	(double*) xmom_explicit_update -> data,
	374	number_of_elements,
	375	(double*) max_speed_array -> data);
[5563]	376
	377
	378	Py_DECREF(neighbours);
	379	Py_DECREF(neighbour_vertices);
	380	Py_DECREF(normals);
	381	Py_DECREF(areas);
	382	Py_DECREF(stage_vertex_values);
	383	Py_DECREF(xmom_vertex_values);
	384	Py_DECREF(bed_vertex_values);
	385	Py_DECREF(stage_boundary_values);
	386	Py_DECREF(xmom_boundary_values);
	387	Py_DECREF(stage_explicit_update);
	388	Py_DECREF(xmom_explicit_update);
	389	Py_DECREF(max_speed_array);
	390
	391
	392
	393
	394	// Return updated flux timestep
	395	return Py_BuildValue("d", timestep);
	396	}
	397
	398
	399
	400
	401	//-------------------------------
	402	// Method table for python module
	403	//-------------------------------
	404
	405	static struct PyMethodDef MethodTable[] = {
	406	{"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
	407	{"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
	408	{NULL, NULL}
	409	};
	410
	411	// Module initialisation
	412	void initcomp_flux_ext(void){
	413	Py_InitModule("comp_flux_ext", MethodTable);
	414	import_array();
	415	}

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

Download in other formats: