Context Navigation

source: anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c @ 7818

Last change on this file since 7818 was 6694, checked in by sudi, 14 years ago
Revised codes for discontinuous bed.
File size: 8.9 KB

Rev	Line
[5827]	1	#include "Python.h"
	2	#include "Numeric/arrayobject.h"
	3	#include "math.h"
	4	#include <stdio.h>
[6694]	5	#include "util_ext.h"
[5827]	6	const double pi = 3.14159265358979;
	7
	8	//Innermost flux function (using w=z+h)
[6694]	9	int _flux_function_wellbalanced(double *q_left,
	10	double *q_right,
	11	double normals,
	12	double g,
	13	double epsilon,
	14	double *edgeflux,
	15	double *max_speed) {
[5827]	16	int i;
	17	double z_left, z_right;
[5832]	18	double flux_left[2], flux_right[2];
[5827]	19	double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
[6694]	20	double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;
[5827]	21	double s_max, s_min, denom;
	22
[6694]	23	w_left = q_left[0];
[5827]	24	uh_left = q_left[1];
	25	uh_left = uh_left*normals;
[6694]	26	z_left = q_left[2];
	27	h_left = q_left[3];
	28	u_left = q_left[4];
[5827]	29
[6694]	30	w_right = q_right[0];
[5827]	31	uh_right = q_right[1];
	32	uh_right = uh_right*normals;
[6694]	33	z_right = q_right[2];
	34	h_right = q_right[3];
	35	u_right = q_right[4];
[5827]	36
	37	if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
	38
	39	z_star = max(z_left, z_right); //equation(7)
	40
[6694]	41	// Compute speeds in x-direction.
	42	h_left_star = max(0.0, w_left-z_star); //equation(8)
	43
[5827]	44	if (h_left_star < epsilon) {
	45	u_left = 0.0; h_left_star = 0.0;
	46	}
	47	else {
	48	u_left = uh_left/h_left_star;
	49	}
[6694]	50
	51	h_right_star = max(0.0, w_right-z_star); //equation(9)
[5827]	52
	53	if (h_right_star < epsilon) {
	54	u_right = 0.0; h_right_star = 0.0;
	55	}
	56	else {
	57	u_right = uh_right/h_right_star;
	58	}
	59
	60	//printf("u_left = %g u_right = %g \n",u_left, u_right);
	61	soundspeed_left = sqrt(g*h_left_star);
	62	soundspeed_right = sqrt(g*h_right_star);
	63
	64	s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
	65	if (s_max < 0.0) s_max = 0.0;
	66
	67	s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
	68	if (s_min > 0.0) s_min = 0.0;
	69
	70
	71	// Flux formulas
	72	flux_left[0] = u_left*h_left_star;
	73	flux_left[1] = u_leftu_lefth_left_star + 0.5gh_left_star*h_left_star;
	74
	75	flux_right[0] = u_right*h_right_star;
	76	flux_right[1] = u_rightu_righth_right_star + 0.5gh_right_star*h_right_star;
	77
	78	// Flux computation
	79	denom = s_max-s_min;
[6694]	80	if (denom < epsilon && z_right > z_left) {
[5827]	81	for (i=0; i<2; i++) edgeflux[i] = 0.0;
[6694]	82	// Maximal wavespeed
	83	edgeflux[1] = edgeflux[1] + (0.5gh_lefth_left - 0.5gh_left_starh_left_star );
	84	edgeflux[1] *= normals;
[5827]	85	*max_speed = 0.0;
[6694]	86	} else if (denom < epsilon){
	87	for (i=0; i<2; i++) edgeflux[i] = 0.0;
	88	// Maximal wavespeed
	89	*max_speed = 0.0;
	90	}else {
[5827]	91	edgeflux[0] = s_maxflux_left[0] - s_minflux_right[0];
	92	edgeflux[0] += s_maxs_min(h_right_star-h_left_star); //s_maxs_min(qr[0]-ql[0]);
	93	edgeflux[0] /= denom;
	94	edgeflux[1] = s_maxflux_left[1] - s_minflux_right[1];
	95	edgeflux[1] += s_maxs_min(flux_right[0]-flux_left[0]); //s_maxs_min(qr[1]-ql[1]);
	96	edgeflux[1] /= denom;
[6694]	97	edgeflux[1] = edgeflux[1] + (0.5gh_lefth_left - 0.5gh_left_starh_left_star );
[5827]	98	edgeflux[1] *= normals;
[6694]	99	// Maximal wavespeed
	100	*max_speed = max(fabs(s_max), fabs(s_min));
[5827]	101	}
[6694]	102
[5827]	103	return 0;
	104	}
	105
	106
	107
	108
	109	// Computational function for flux computation
	110	double _compute_fluxes_ext_wellbalanced(double timestep,
	111	double epsilon,
	112	double g,
[6694]	113
[5827]	114	long* neighbours,
	115	long* neighbour_vertices,
	116	double* normals,
	117	double* areas,
[6694]	118
[5827]	119	double* stage_edge_values,
	120	double* xmom_edge_values,
	121	double* bed_edge_values,
[6694]	122	double* height_edge_values,
	123	double* vel_edge_values,
	124
[5827]	125	double* stage_boundary_values,
	126	double* xmom_boundary_values,
[6694]	127	double* bed_boundary_values,
	128	double* height_boundary_values,
	129	double* vel_boundary_values,
	130
[5827]	131	double* stage_explicit_update,
	132	double* xmom_explicit_update,
[6694]	133	double* bed_explicit_update,
	134	double* height_explicit_update,
	135	double* vel_explicit_update,
	136
[5827]	137	int number_of_elements,
	138	double* max_speed_array) {
	139
[6694]	140	double flux[2], ql[5], qr[5], edgeflux[2];
[5832]	141	double max_speed, normal;
[5827]	142	int k, i, ki, n, m, nm=0;
	143
	144	for (k=0; k<number_of_elements; k++) {
	145	flux[0] = 0.0;
	146	flux[1] = 0.0;
	147
	148	for (i=0; i<2; i++) {
	149	ki = k*2+i;
	150
	151	ql[0] = stage_edge_values[ki];
	152	ql[1] = xmom_edge_values[ki];
	153	ql[2] = bed_edge_values[ki];
[6694]	154	ql[3] = height_edge_values[ki];
	155	ql[4] = vel_edge_values[ki];
[5827]	156
	157	n = neighbours[ki];
[6694]	158	printf("neighbours[ki] = %li \n",neighbours[ki]);
[5827]	159	if (n<0) {
	160	m = -n-1;
	161	qr[0] = stage_boundary_values[m];
	162	qr[1] = xmom_boundary_values[m];
[6694]	163	qr[2] = ql[2];
	164	qr[3] = height_edge_values[m];
	165	qr[4] = vel_edge_values[m];
[5827]	166
	167	} else {
	168	m = neighbour_vertices[ki];
	169	nm = n*2+m;
	170	qr[0] = stage_edge_values[nm];
	171	qr[1] = xmom_edge_values[nm];
[6694]	172	qr[2] = bed_edge_values[nm];
	173	qr[3] = height_edge_values[nm];
	174	qr[4] = vel_edge_values[nm];
[5827]	175	}
	176
	177	normal = normals[ki];
	178	_flux_function_wellbalanced(ql, qr, normal, g, epsilon, edgeflux, &max_speed);
	179	flux[0] -= edgeflux[0];
	180	flux[1] -= edgeflux[1];
	181
	182	// Update timestep based on edge i and possibly neighbour n
	183	if (max_speed > epsilon) {
	184	// Original CFL calculation
[6694]	185	timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed.
[5827]	186	if (n>=0) {
[6694]	187	timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed.
[5827]	188	}
	189	}
	190	} // End edge i (and neighbour n)
	191	flux[0] /= areas[k];
	192	stage_explicit_update[k] = flux[0];
	193	flux[1] /= areas[k];
	194	xmom_explicit_update[k] = flux[1];
	195
	196	//Keep track of maximal speeds
	197	max_speed_array[k]=max_speed;
	198	}
[6694]	199
[5827]	200	return timestep;
	201	}
	202
	203
	204
	205
	206
	207
	208	//=========================================================================
	209	// Python Glue
	210	//=========================================================================
	211	PyObject compute_fluxes_ext_wellbalanced(PyObject self, PyObject *args) {
	212
	213	PyArrayObject
	214	*neighbours,
	215	*neighbour_vertices,
	216	*normals,
	217	*areas,
[6694]	218
[5827]	219	*stage_edge_values,
	220	*xmom_edge_values,
	221	*bed_edge_values,
[6694]	222	*height_edge_values,
	223	*vel_edge_values,
	224
[5827]	225	*stage_boundary_values,
	226	*xmom_boundary_values,
[6694]	227	*bed_boundary_values,
	228	*height_boundary_values,
	229	*vel_boundary_values,
	230
[5827]	231	*stage_explicit_update,
	232	*xmom_explicit_update,
[6694]	233	*bed_explicit_update,
	234	*height_explicit_update,
	235	*vel_explicit_update,
	236
[5827]	237	*max_speed_array;
	238
	239	double timestep, epsilon, g;
	240	int number_of_elements;
	241
	242	// Convert Python arguments to C
[6694]	243	if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOOiO",
[5827]	244	&timestep,
	245	&epsilon,
	246	&g,
[6694]	247
[5827]	248	&neighbours,
	249	&neighbour_vertices,
	250	&normals,
	251	&areas,
[6694]	252
[5827]	253	&stage_edge_values,
	254	&xmom_edge_values,
	255	&bed_edge_values,
[6694]	256	&height_edge_values,
	257	&vel_edge_values,
	258
[5827]	259	&stage_boundary_values,
	260	&xmom_boundary_values,
[6694]	261	&bed_boundary_values,
	262	&height_boundary_values,
	263	&vel_boundary_values,
	264
[5827]	265	&stage_explicit_update,
	266	&xmom_explicit_update,
[6694]	267	&bed_explicit_update,
	268	&height_explicit_update,
	269	&vel_explicit_update,
	270
[5827]	271	&number_of_elements,
	272	&max_speed_array)) {
[6694]	273	PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext_wellbalanced.c: compute_fluxes_ext_wellbalanced could not parse input");
[5827]	274	return NULL;
	275	}
	276
	277
	278	// Call underlying flux computation routine and update
	279	// the explicit update arrays
	280	timestep = _compute_fluxes_ext_wellbalanced(timestep,
	281	epsilon,
	282	g,
[6694]	283
[5827]	284	(long*) neighbours -> data,
	285	(long*) neighbour_vertices -> data,
	286	(double*) normals -> data,
	287	(double*) areas -> data,
[6694]	288
[5827]	289	(double*) stage_edge_values -> data,
	290	(double*) xmom_edge_values -> data,
	291	(double*) bed_edge_values -> data,
[6694]	292	(double*) height_edge_values -> data,
	293	(double*) vel_edge_values -> data,
	294
[5827]	295	(double*) stage_boundary_values -> data,
	296	(double*) xmom_boundary_values -> data,
[6694]	297	(double*) bed_boundary_values -> data,
	298	(double*) height_boundary_values -> data,
	299	(double*) vel_boundary_values -> data,
	300
[5827]	301	(double*) stage_explicit_update -> data,
	302	(double*) xmom_explicit_update -> data,
[6694]	303	(double*) bed_explicit_update -> data,
	304	(double*) height_explicit_update -> data,
	305	(double*) vel_explicit_update -> data,
	306
[5827]	307	number_of_elements,
	308	(double*) max_speed_array -> data);
	309
	310	// Return updated flux timestep
	311	return Py_BuildValue("d", timestep);
	312	}
	313
	314
	315
	316	//-------------------------------
	317	// Method table for python module
	318	//-------------------------------
	319
	320	static struct PyMethodDef MethodTable[] = {
	321	{"compute_fluxes_ext_wellbalanced", compute_fluxes_ext_wellbalanced, METH_VARARGS, "Print out"},
	322	{NULL, NULL}
	323	};
	324
	325	// Module initialisation
	326	void initcomp_flux_ext_wellbalanced(void){
	327	Py_InitModule("comp_flux_ext_wellbalanced", MethodTable);
	328	import_array();
	329	}

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

Download in other formats: