Changeset 4681


Ignore:
Timestamp:
Aug 23, 2007, 4:37:38 PM (17 years ago)
Author:
ole
Message:

Cleaned up flux computation code and profiled.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4677 r4681  
    180180// Computational function for flux computation (using stage w=z+h)
    181181int 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) {
    186186
    187187  /*Compute fluxes between volumes for the shallow water wave equation
     
    226226  h_left = w_left-z;
    227227  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);
    229229
    230230  w_right = q_right_copy[0];
    231231  h_right = w_right-z;
    232232  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); 
    234234
    235235  //Momentum in y-direction
     
    237237  vh_right = q_right_copy[2];
    238238
    239   // Limit momentum if necessary 
    240   v_left =_compute_speed(&vh_left, &h_left, epsilon, h0);
    241   v_right =_compute_speed(&vh_right, &h_right, epsilon, h0);
     239  // Limit y-momentum if necessary 
     240  v_left = _compute_speed(&vh_left, &h_left, epsilon, h0);
     241  v_right = _compute_speed(&vh_right, &h_right, epsilon, h0);
    242242
    243243  //Maximal and minimal wave speeds
     
    280280    _rotate(edgeflux, n1, -n2);
    281281  }
     282 
    282283  return 0;
    283284}
     
    818819    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
    819820    original triangle has gradient (a,b).
    820     3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
     821    3) Provisional vertex jumps dqv[0,1,2] are computed and these are then limited by calling the functions
    821822    find_qmin_and_qmax and limit_gradient
    822823
     
    13061307    *max_speed_array; //Keeps track of max speeds for each triangle
    13071308
    1308   //Local variables
    1309   double timestep, max_speed, epsilon, g, H0;
     1309  // Local variables
     1310  double timestep, max_speed, epsilon, g, H0, length;
    13101311  double normal[2], ql[3], qr[3], zl, zr;
    1311   double edgeflux[3]; //Work arrays for summing up fluxes
     1312  double edgeflux[3]; // Work array for summing up fluxes
    13121313
    13131314  int number_of_elements, k, i, m, n;
    1314   int ki, nm=0, ki2; //Index shorthands
    1315   static long call=1;
     1315  int ki, nm=0, ki2; // Index shorthands
     1316  static long call=1; // Static local variable flagging already computed flux
    13161317
    13171318
     
    13421343    return NULL;
    13431344  }
     1345 
    13441346  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
    13481352  for (k=0; k<number_of_elements; k++) {
    13491353    ((double *) stage_explicit_update -> data)[k]=0.0;
     
    13521356  }
    13531357 
    1354   //Loop through neighbours and compute edge flux for each
     1358  // For all triangles
    13551359  for (k=0; k<number_of_elements; k++) {
     1360 
     1361    // Loop through neighbours and compute edge flux for each 
    13561362    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
    13591367        continue;
     1368       
    13601369      ql[0] = ((double *) stage_edge_values -> data)[ki];
    13611370      ql[1] = ((double *) xmom_edge_values -> data)[ki];
     
    13631372      zl =    ((double *) bed_edge_values -> data)[ki];
    13641373
    1365       //Quantities at neighbour on nearest face
     1374      // Quantities at neighbour on nearest face
    13661375      n = ((long *) neighbours -> data)[ki];
    13671376      if (n < 0) {
    1368         m = -n-1; //Convert negative flag to index
     1377        m = -n-1; // Convert negative flag to index
     1378       
    13691379        qr[0] = ((double *) stage_boundary_values -> data)[m];
    13701380        qr[1] = ((double *) xmom_boundary_values -> data)[m];
     
    13731383      } else {
    13741384        m = ((long *) neighbour_edges -> data)[ki];
    1375         nm = n*3+m;
     1385        nm = n*3+m; // Linear index (triangle n, edge m)
     1386       
    13761387        qr[0] = ((double *) stage_edge_values -> data)[nm];
    13771388        qr[1] = ((double *) xmom_edge_values -> data)[nm];
     
    13791390        zr =    ((double *) bed_edge_values -> data)[nm];
    13801391      }
    1381       // Outward pointing normal vector
    1382       // normal = domain.normals[k, 2*i:2*i+2]
     1392     
     1393      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
    13831394      ki2 = 2*ki; //k*6 + i*2
    13841395      normal[0] = ((double *) normals -> data)[ki2];
    13851396      normal[1] = ((double *) normals -> data)[ki2+1];
    1386       //Edge flux computation
     1397     
     1398      // Edge flux computation (triangle k, edge i)
    13871399      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
    13971421      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|,|u-a|)
    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
    14181439   
    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
    14211443    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
    14221444    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
     
    14241446   
    14251447   
    1426     //Keep track of maximal speeds
     1448    // Keep track of maximal speeds
    14271449    ((double *) max_speed_array -> data)[k] = max_speed;   
    14281450   
    1429   } //end for k
     1451  } // End triangle k
     1452 
    14301453  return Py_BuildValue("d", timestep);
    14311454}
Note: See TracChangeset for help on using the changeset viewer.