Changeset 4681
 Timestamp:
 Aug 23, 2007, 4:37:38 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r4677 r4681 180 180 // Computational function for flux computation (using stage w=z+h) 181 181 int flux_function_central(double *q_left, double *q_right, 182 double z_left, double z_right,183 double n1, double n2,184 double epsilon, double H0, double g,185 double *edgeflux, double *max_speed) {182 double z_left, double z_right, 183 double n1, double n2, 184 double epsilon, double H0, double g, 185 double *edgeflux, double *max_speed) { 186 186 187 187 /*Compute fluxes between volumes for the shallow water wave equation … … 226 226 h_left = w_leftz; 227 227 uh_left = q_left_copy[1]; 228 u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);228 u_left = _compute_speed(&uh_left, &h_left, epsilon, h0); 229 229 230 230 w_right = q_right_copy[0]; 231 231 h_right = w_rightz; 232 232 uh_right = q_right_copy[1]; 233 u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);233 u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); 234 234 235 235 //Momentum in ydirection … … 237 237 vh_right = q_right_copy[2]; 238 238 239 // Limit momentum if necessary240 v_left = _compute_speed(&vh_left, &h_left, epsilon, h0);241 v_right = _compute_speed(&vh_right, &h_right, epsilon, h0);239 // Limit ymomentum if necessary 240 v_left = _compute_speed(&vh_left, &h_left, epsilon, h0); 241 v_right = _compute_speed(&vh_right, &h_right, epsilon, h0); 242 242 243 243 //Maximal and minimal wave speeds … … 280 280 _rotate(edgeflux, n1, n2); 281 281 } 282 282 283 return 0; 283 284 } … … 818 819 of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the 819 820 original triangle has gradient (a,b). 820 3) Provisional vertex ju nmps dqv[0,1,2] are computed and these are then limited by calling the functions821 3) Provisional vertex jumps dqv[0,1,2] are computed and these are then limited by calling the functions 821 822 find_qmin_and_qmax and limit_gradient 822 823 … … 1306 1307 *max_speed_array; //Keeps track of max speeds for each triangle 1307 1308 1308 // Local variables1309 double timestep, max_speed, epsilon, g, H0 ;1309 // Local variables 1310 double timestep, max_speed, epsilon, g, H0, length; 1310 1311 double normal[2], ql[3], qr[3], zl, zr; 1311 double edgeflux[3]; // Work arraysfor summing up fluxes1312 double edgeflux[3]; // Work array for summing up fluxes 1312 1313 1313 1314 int number_of_elements, k, i, m, n; 1314 int ki, nm=0, ki2; // Index shorthands1315 static long call=1; 1315 int ki, nm=0, ki2; // Index shorthands 1316 static long call=1; // Static local variable flagging already computed flux 1316 1317 1317 1318 … … 1342 1343 return NULL; 1343 1344 } 1345 1344 1346 number_of_elements = stage_edge_values > dimensions[0]; 1345 call++;//a static local variable to which already_computed_flux is compared 1346 //set explicit_update to zero for all conserved_quantities. 1347 //This assumes compute_fluxes called before forcing terms 1347 1348 call++; // Flag 'id' of flux calculation for this timestep 1349 1350 // Set explicit_update to zero for all conserved_quantities. 1351 // This assumes compute_fluxes called before forcing terms 1348 1352 for (k=0; k<number_of_elements; k++) { 1349 1353 ((double *) stage_explicit_update > data)[k]=0.0; … … 1352 1356 } 1353 1357 1354 // Loop through neighbours and compute edge flux for each1358 // For all triangles 1355 1359 for (k=0; k<number_of_elements; k++) { 1360 1361 // Loop through neighbours and compute edge flux for each 1356 1362 for (i=0; i<3; i++) { 1357 ki = k*3+i; 1358 if (((long *) already_computed_flux>data)[ki]==call)//we've already computed the flux across this edge 1363 ki = k*3+i; // Linear index (triangle k, edge i) 1364 1365 if (((long *) already_computed_flux>data)[ki] == call) 1366 // We've already computed the flux across this edge 1359 1367 continue; 1368 1360 1369 ql[0] = ((double *) stage_edge_values > data)[ki]; 1361 1370 ql[1] = ((double *) xmom_edge_values > data)[ki]; … … 1363 1372 zl = ((double *) bed_edge_values > data)[ki]; 1364 1373 1365 // Quantities at neighbour on nearest face1374 // Quantities at neighbour on nearest face 1366 1375 n = ((long *) neighbours > data)[ki]; 1367 1376 if (n < 0) { 1368 m = n1; //Convert negative flag to index 1377 m = n1; // Convert negative flag to index 1378 1369 1379 qr[0] = ((double *) stage_boundary_values > data)[m]; 1370 1380 qr[1] = ((double *) xmom_boundary_values > data)[m]; … … 1373 1383 } else { 1374 1384 m = ((long *) neighbour_edges > data)[ki]; 1375 nm = n*3+m; 1385 nm = n*3+m; // Linear index (triangle n, edge m) 1386 1376 1387 qr[0] = ((double *) stage_edge_values > data)[nm]; 1377 1388 qr[1] = ((double *) xmom_edge_values > data)[nm]; … … 1379 1390 zr = ((double *) bed_edge_values > data)[nm]; 1380 1391 } 1381 // Outward pointing normal vector1382 // normal = domain.normals[k, 2*i:2*i+2]1392 1393 // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2]) 1383 1394 ki2 = 2*ki; //k*6 + i*2 1384 1395 normal[0] = ((double *) normals > data)[ki2]; 1385 1396 normal[1] = ((double *) normals > data)[ki2+1]; 1386 //Edge flux computation 1397 1398 // Edge flux computation (triangle k, edge i) 1387 1399 flux_function_central(ql, qr, zl, zr, 1388 normal[0], normal[1], 1389 epsilon, H0, g, 1390 edgeflux, &max_speed); 1391 //update triangle k 1392 ((long *) already_computed_flux>data)[ki]=call; 1393 ((double *) stage_explicit_update > data)[k] = edgeflux[0]*((double *) edgelengths > data)[ki]; 1394 ((double *) xmom_explicit_update > data)[k] = edgeflux[1]*((double *) edgelengths > data)[ki]; 1395 ((double *) ymom_explicit_update > data)[k] = edgeflux[2]*((double *) edgelengths > data)[ki]; 1396 //update the neighbour n 1400 normal[0], normal[1], 1401 epsilon, H0, g, 1402 edgeflux, &max_speed); 1403 1404 1405 // Multiply edgeflux by edgelength 1406 length = ((double *) edgelengths > data)[ki]; 1407 edgeflux[0] *= length; 1408 edgeflux[1] *= length; 1409 edgeflux[2] *= length; 1410 1411 1412 // Update triangle k with flux from edge i 1413 ((double *) stage_explicit_update > data)[k] = edgeflux[0]; 1414 ((double *) xmom_explicit_update > data)[k] = edgeflux[1]; 1415 ((double *) ymom_explicit_update > data)[k] = edgeflux[2]; 1416 1417 ((long *) already_computed_flux > data)[ki] = call; // #k Done 1418 1419 1420 // Update neighbour n with same flux but reversed sign 1397 1421 if (n>=0){ 1398 ((long *) already_computed_flux>data)[nm]=call; 1399 ((double *) stage_explicit_update > data)[n] += edgeflux[0]*((double *) edgelengths > data)[nm]; 1400 ((double *) xmom_explicit_update > data)[n] += edgeflux[1]*((double *) edgelengths > data)[nm]; 1401 ((double *) ymom_explicit_update > data)[n] += edgeflux[2]*((double *) edgelengths > data)[nm]; 1402 } 1403 ///for (j=0; j<3; j++) { 1404 ///flux[j] = edgeflux[j]*((double *) edgelengths > data)[ki]; 1405 ///} 1406 //Update timestep 1407 //timestep = min(timestep, domain.radii[k]/max_speed) 1408 //FIXME: SR Add parameter for CFL condition 1409 if ( ((long *) tri_full_flag > data)[k] == 1) { 1410 if (max_speed > epsilon) { 1411 timestep = min(timestep, ((double *) radii > data)[k]/max_speed); 1412 //maxspeed in flux_function is calculated as max(u+a,ua) 1413 if (n>=0) 1414 timestep = min(timestep, ((double *) radii > data)[n]/max_speed); 1415 } 1416 } 1417 } // end for i 1422 ((double *) stage_explicit_update > data)[n] += edgeflux[0]; 1423 ((double *) xmom_explicit_update > data)[n] += edgeflux[1]; 1424 ((double *) ymom_explicit_update > data)[n] += edgeflux[2]; 1425 1426 ((long *) already_computed_flux > data)[nm] = call; // #n Done 1427 } 1428 1429 1430 // Update timestep based on edge i and possibly neighbour n 1431 if ( ((long *) tri_full_flag > data)[k] == 1) { 1432 if (max_speed > epsilon) { 1433 timestep = min(timestep, ((double *) radii > data)[k]/max_speed); 1434 if (n>=0) 1435 timestep = min(timestep, ((double *) radii > data)[n]/max_speed); 1436 } 1437 } 1438 } // End edge i 1418 1439 1419 //Normalise by area and store for when all conserved 1420 //quantities get updated 1440 1441 // Normalise triangle k by area and store for when all conserved 1442 // quantities get updated 1421 1443 ((double *) stage_explicit_update > data)[k] /= ((double *) areas > data)[k]; 1422 1444 ((double *) xmom_explicit_update > data)[k] /= ((double *) areas > data)[k]; … … 1424 1446 1425 1447 1426 // Keep track of maximal speeds1448 // Keep track of maximal speeds 1427 1449 ((double *) max_speed_array > data)[k] = max_speed; 1428 1450 1429 } //end for k 1451 } // End triangle k 1452 1430 1453 return Py_BuildValue("d", timestep); 1431 1454 }
Note: See TracChangeset
for help on using the changeset viewer.