- Timestamp:
- Apr 1, 2009, 8:41:37 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c
r5832 r6694 3 3 #include "math.h" 4 4 #include <stdio.h> 5 #include "util_ext.h" 5 6 const double pi = 3.14159265358979; 6 7 7 8 // Shared code snippets9 #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 24 25 26 8 //Innermost flux function (using w=z+h) 27 int _flux_function_wellbalanced(double *q_left, double *q_right, 28 double normals, double g, double epsilon, 29 double *edgeflux, double *max_speed) { 30 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) { 31 16 int i; 32 17 double z_left, z_right; 33 18 double flux_left[2], flux_right[2]; 34 19 double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left; 35 double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right; //h_right,20 double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right; 36 21 double s_max, s_min, denom; 37 22 38 w_left = q_left[0];23 w_left = q_left[0]; 39 24 uh_left = q_left[1]; 40 25 uh_left = uh_left*normals; 41 z_left = q_left[2]; 42 43 w_right = q_right[0]; 26 z_left = q_left[2]; 27 h_left = q_left[3]; 28 u_left = q_left[4]; 29 30 w_right = q_right[0]; 44 31 uh_right = q_right[1]; 45 32 uh_right = uh_right*normals; 46 z_right = q_right[2]; 33 z_right = q_right[2]; 34 h_right = q_right[3]; 35 u_right = q_right[4]; 47 36 48 37 if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right); 49 38 50 //z = (z_left+z_right)/2.0;51 39 z_star = max(z_left, z_right); //equation(7) 52 40 53 // Compute speeds in x-direction 54 h_left = w_left-z_left; //This is used for computing the edgeflux[1]. 55 h_left_star = max(0, w_left-z_star); //(8) 56 41 // Compute speeds in x-direction. 42 h_left_star = max(0.0, w_left-z_star); //equation(8) 43 57 44 if (h_left_star < epsilon) { 58 45 u_left = 0.0; h_left_star = 0.0; … … 61 48 u_left = uh_left/h_left_star; 62 49 } 63 64 h_right = w_right-z_right; 65 h_right_star = max(0, w_right-z_star); //(9) 50 51 h_right_star = max(0.0, w_right-z_star); //equation(9) 66 52 67 53 if (h_right_star < epsilon) { … … 92 78 // Flux computation 93 79 denom = s_max-s_min; 94 if (denom < epsilon ) {80 if (denom < epsilon && z_right > z_left) { 95 81 for (i=0; i<2; i++) edgeflux[i] = 0.0; 82 // Maximal wavespeed 83 edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star ); 84 edgeflux[1] *= normals; 96 85 *max_speed = 0.0; 97 } else { 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 { 98 91 edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0]; 99 92 edgeflux[0] += s_max*s_min*(h_right_star-h_left_star); //s_max*s_min*(qr[0]-ql[0]); … … 102 95 edgeflux[1] += s_max*s_min*(flux_right[0]-flux_left[0]); //s_max*s_min*(qr[1]-ql[1]); 103 96 edgeflux[1] /= denom; 97 edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star ); 104 98 edgeflux[1] *= normals; 105 106 //edgeflux[1] = edgeflux[1] + 0.5*g*0.5*(h_left+h_right)*0.5*(h_left+h_right) - 0.5*g*h_left_star*h_left_star; 107 //edgeflux[1] = edgeflux[1] + 0.5*g*h_right*h_right - 0.5*g*h_right_star*h_right_star; 108 edgeflux[1] = edgeflux[1] + 0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star; //what about the h_right and h_right_star?????????????????????? 109 // Maximal wavespeed 110 *max_speed = max(fabs(s_max), fabs(s_min)); 99 // Maximal wavespeed 100 *max_speed = max(fabs(s_max), fabs(s_min)); 111 101 } 112 //printf("max_speed = %g \n",*max_speed);102 113 103 return 0; 114 104 } … … 121 111 double epsilon, 122 112 double g, 113 123 114 long* neighbours, 124 115 long* neighbour_vertices, 125 116 double* normals, 126 117 double* areas, 118 127 119 double* stage_edge_values, 128 120 double* xmom_edge_values, 129 121 double* bed_edge_values, 122 double* height_edge_values, 123 double* vel_edge_values, 124 130 125 double* stage_boundary_values, 131 126 double* xmom_boundary_values, 127 double* bed_boundary_values, 128 double* height_boundary_values, 129 double* vel_boundary_values, 130 132 131 double* stage_explicit_update, 133 132 double* xmom_explicit_update, 133 double* bed_explicit_update, 134 double* height_explicit_update, 135 double* vel_explicit_update, 136 134 137 int number_of_elements, 135 138 double* max_speed_array) { 136 139 137 double flux[2], ql[ 3], qr[3], edgeflux[2];140 double flux[2], ql[5], qr[5], edgeflux[2]; 138 141 double max_speed, normal; 139 142 int k, i, ki, n, m, nm=0; … … 149 152 ql[1] = xmom_edge_values[ki]; 150 153 ql[2] = bed_edge_values[ki]; 151 //ql[3] = velocity_edge_values[ki];152 //ql[4] = height_edge_values[ki];154 ql[3] = height_edge_values[ki]; 155 ql[4] = vel_edge_values[ki]; 153 156 154 157 n = neighbours[ki]; 158 printf("neighbours[ki] = %li \n",neighbours[ki]); 155 159 if (n<0) { 156 160 m = -n-1; 157 161 qr[0] = stage_boundary_values[m]; 158 162 qr[1] = xmom_boundary_values[m]; 159 qr[2] = ql[2]; 163 qr[2] = ql[2]; 164 qr[3] = height_edge_values[m]; 165 qr[4] = vel_edge_values[m]; 160 166 161 167 } else { … … 164 170 qr[0] = stage_edge_values[nm]; 165 171 qr[1] = xmom_edge_values[nm]; 166 qr[2] = bed_edge_values[nm]; 172 qr[2] = bed_edge_values[nm]; 173 qr[3] = height_edge_values[nm]; 174 qr[4] = vel_edge_values[nm]; 167 175 } 168 176 … … 175 183 if (max_speed > epsilon) { 176 184 // Original CFL calculation 177 178 timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ????????????????????????????????????????????? 185 timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. 179 186 if (n>=0) { 180 timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????187 timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. 181 188 } 182 189 } … … 190 197 max_speed_array[k]=max_speed; 191 198 } 192 //printf("timestep = %f \n",timestep);199 193 200 return timestep; 194 201 } … … 209 216 *normals, 210 217 *areas, 218 211 219 *stage_edge_values, 212 220 *xmom_edge_values, 213 221 *bed_edge_values, 222 *height_edge_values, 223 *vel_edge_values, 224 214 225 *stage_boundary_values, 215 226 *xmom_boundary_values, 227 *bed_boundary_values, 228 *height_boundary_values, 229 *vel_boundary_values, 230 216 231 *stage_explicit_update, 217 232 *xmom_explicit_update, 233 *bed_explicit_update, 234 *height_explicit_update, 235 *vel_explicit_update, 236 218 237 *max_speed_array; 219 238 … … 222 241 223 242 // Convert Python arguments to C 224 if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOO iO",243 if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOOiO", 225 244 ×tep, 226 245 &epsilon, 227 246 &g, 247 228 248 &neighbours, 229 249 &neighbour_vertices, 230 250 &normals, 231 251 &areas, 252 232 253 &stage_edge_values, 233 254 &xmom_edge_values, 234 255 &bed_edge_values, 256 &height_edge_values, 257 &vel_edge_values, 258 235 259 &stage_boundary_values, 236 260 &xmom_boundary_values, 261 &bed_boundary_values, 262 &height_boundary_values, 263 &vel_boundary_values, 264 237 265 &stage_explicit_update, 238 266 &xmom_explicit_update, 267 &bed_explicit_update, 268 &height_explicit_update, 269 &vel_explicit_update, 270 239 271 &number_of_elements, 240 272 &max_speed_array)) { 241 PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext .c: compute_fluxes_extcould not parse input");273 PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext_wellbalanced.c: compute_fluxes_ext_wellbalanced could not parse input"); 242 274 return NULL; 243 275 } … … 249 281 epsilon, 250 282 g, 283 251 284 (long*) neighbours -> data, 252 285 (long*) neighbour_vertices -> data, 253 286 (double*) normals -> data, 254 287 (double*) areas -> data, 288 255 289 (double*) stage_edge_values -> data, 256 290 (double*) xmom_edge_values -> data, 257 291 (double*) bed_edge_values -> data, 292 (double*) height_edge_values -> data, 293 (double*) vel_edge_values -> data, 294 258 295 (double*) stage_boundary_values -> data, 259 296 (double*) xmom_boundary_values -> data, 297 (double*) bed_boundary_values -> data, 298 (double*) height_boundary_values -> data, 299 (double*) vel_boundary_values -> data, 300 260 301 (double*) stage_explicit_update -> data, 261 302 (double*) xmom_explicit_update -> data, 303 (double*) bed_explicit_update -> data, 304 (double*) height_explicit_update -> data, 305 (double*) vel_explicit_update -> data, 306 262 307 number_of_elements, 263 308 (double*) max_speed_array -> data); … … 266 311 return Py_BuildValue("d", timestep); 267 312 } 268 269 270 313 271 314 … … 285 328 import_array(); 286 329 } 287 288
Note: See TracChangeset
for help on using the changeset viewer.