Ignore:
Timestamp:
Jun 9, 2005, 3:37:55 PM (19 years ago)
Author:
matthew
Message:

To avoid computing the flux function at each internal edge twice, we now call the C implementation of compute_fluxes with the extra integer vector argument already_computed_flux. This is an array of dimension (N,3) in python (or 3N in C) and thus each entry corresponds to a triangle-edge pair. If triangle t has neighbour n across its i-th edge, and this edge is the m-th edge of triangle n, then we update this edge's flux contribution to explicit_update for both triangle n and triangle t. Both corresponding entries of already_computed_flux will be incremented so that, when we come to edge m of triangle n, we know that the update across this edge has already been done.
The array already_computed_flux is defined and initialised to zero in domain.py.

For run_profile.py, with N=128, the version of compute_fluxes in which each edge was computed twice consumed 13.64 seconds (averaged over 3 runs). By keeping track of which edges have been computed, compute_fluxes now consumes 8.68 seconds.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

    r1507 r1508  
    848848}
    849849
    850 
    851850PyObject *compute_fluxes(PyObject *self, PyObject *args) {
    852851  /*Compute all fluxes and the timestep suitable for all volumes
    … …  
    883882                                     Stage.explicit_update,
    884883                                     Xmom.explicit_update,
    885                                      Ymom.explicit_update)
     884                                     Ymom.explicit_update,
     885                                     already_computed_flux)
    886886
    887887
    … …  
    905905    *stage_explicit_update,
    906906    *xmom_explicit_update,
    907     *ymom_explicit_update;
     907    *ymom_explicit_update,
     908    *already_computed_flux;//tracks whether the flux across an edge has already been computed
    908909
    909910
    … …  
    911912  double timestep, max_speed, epsilon, g;
    912913  double normal[2], ql[3], qr[3], zl, zr;
    913   double flux[3], edgeflux[3]; //Work arrays for summing up fluxes
    914 
    915   int number_of_elements, k, i, j, m, n;
     914  double edgeflux[3]; //Work arrays for summing up fluxes
     915
     916  int number_of_elements, k, i, m, n;
    916917  int ki, nm, ki2; //Index shorthands
     918  static long call=1;
    917919
    918920
    919921  // Convert Python arguments to C
    920   if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOO",
     922  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOO",
    921923                        &timestep,
    922924                        &epsilon,
    … …  
    935937                        &stage_explicit_update,
    936938                        &xmom_explicit_update,
    937                         &ymom_explicit_update)) {
     939                        &ymom_explicit_update,
     940                        &already_computed_flux)) {
    938941    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    939942    return NULL;
    940943  }
    941 
    942944  number_of_elements = stage_edge_values -> dimensions[0];
    943 
    944 
     945  call++;//a static local variable to which already_computed_flux is compared
     946  //set explicit_update to zero for all conserved_quantities.
     947  //This assumes compute_fluxes called before forcing terms
    945948  for (k=0; k<number_of_elements; k++) {
    946 
    947     //Reset work array
    948     for (j=0; j<3; j++) flux[j] = 0.0;
    949 
    950     //Loop through neighbours and compute edge flux for each
     949    ((double *) stage_explicit_update -> data)[k]=0.0;
     950    ((double *) xmom_explicit_update -> data)[k]=0.0;
     951    ((double *) ymom_explicit_update -> data)[k]=0.0;
     952  }
     953  //Loop through neighbours and compute edge flux for each
     954  for (k=0; k<number_of_elements; k++) {
    951955    for (i=0; i<3; i++) {
    952956      ki = k*3+i;
     957      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
     958        continue;
    953959      ql[0] = ((double *) stage_edge_values -> data)[ki];
    954960      ql[1] = ((double *) xmom_edge_values -> data)[ki];
    955961      ql[2] = ((double *) ymom_edge_values -> data)[ki];
    956962      zl =    ((double *) bed_edge_values -> data)[ki];
    957 
     963     
    958964      //Quantities at neighbour on nearest face
    959965      n = ((long *) neighbours -> data)[ki];
    … …  
    966972      } else {
    967973        m = ((long *) neighbour_edges -> data)[ki];
    968 
    969974        nm = n*3+m;
    970975        qr[0] = ((double *) stage_edge_values -> data)[nm];
    … …  
    973978        zr =    ((double *) bed_edge_values -> data)[nm];
    974979      }
    975 
    976       //printf("%d %d [%d] (%d, %d): %.2f %.2f %.2f\n", k, i, ki, n, m,
    977       //             ql[0], ql[1], ql[2]);
    978 
    979 
    980980      // Outward pointing normal vector
    981981      // normal = domain.normals[k, 2*i:2*i+2]
    … …  
    983983      normal[0] = ((double *) normals -> data)[ki2];
    984984      normal[1] = ((double *) normals -> data)[ki2+1];
    985 
    986985      //Edge flux computation
    987986      flux_function(ql, qr, zl, zr,
    … …  
    989988                    epsilon, g,
    990989                    edgeflux, &max_speed);
    991 
    992 
    993       //flux -= edgeflux * edgelengths[k,i]
    994       for (j=0; j<3; j++) {
    995         flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
     990      //update triangle k
     991      ((long *) already_computed_flux->data)[ki]=call;
     992      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
     993      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
     994      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
     995      //update the neighbour n
     996      if (n>=0){
     997        ((long *) already_computed_flux->data)[nm]=call;
     998        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
     999        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
     1000        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
    9961001      }
    997 
    998       //Update timestep
    999       //timestep = min(timestep, domain.radii[k]/max_speed)
    1000       //FIXME: SR Add parameter for CFL condition
    1001       if (max_speed > epsilon) {
    1002         timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
    1003       }
     1002      ///for (j=0; j<3; j++) {
     1003        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
     1004        ///}
     1005        //Update timestep
     1006        //timestep = min(timestep, domain.radii[k]/max_speed)
     1007        //FIXME: SR Add parameter for CFL condition
     1008        if (max_speed > epsilon) {
     1009          timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
     1010          //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
     1011          if (n>=0)
     1012            timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
     1013        }
    10041014    } // end for i
    1005 
    10061015    //Normalise by area and store for when all conserved
    10071016    //quantities get updated
    1008     // flux /= areas[k]
    1009     for (j=0; j<3; j++) {
    1010       flux[j] /= ((double *) areas -> data)[k];
    1011     }
    1012 
    1013     ((double *) stage_explicit_update -> data)[k] = flux[0];
    1014     ((double *) xmom_explicit_update -> data)[k] = flux[1];
    1015     ((double *) ymom_explicit_update -> data)[k] = flux[2];
    1016 
     1017    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
     1018    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
     1019    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];   
    10171020  } //end for k
    1018 
    10191021  return Py_BuildValue("d", timestep);
    10201022}
    1021 
    1022 
    10231023
    10241024PyObject *protect(PyObject *self, PyObject *args) {
Note: See TracChangeset for help on using the changeset viewer.