Changeset 4710


Ignore:
Timestamp:
Sep 6, 2007, 4:12:42 PM (18 years ago)
Author:
ole
Message:

Implemented dry-cell exclusion for the linear reconstruction step.
Time for the okushiri performance example reduced the running time (on my GA windows box)
of that step from 15.3s to 11.6s - an almost 25% improvement.
With a total running time of this example of about 33s before the optimisation this translates to just over 10% total improvement of the running time of the evolve loop.

With more dry land, this optimisation will be more important.

Location:
anuga_core/source/anuga/shallow_water
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r4701 r4710  
    9999from anuga.config import minimum_storable_height
    100100from anuga.config import minimum_allowed_height, maximum_allowed_speed
    101 from anuga.config import g, beta_h, beta_w, beta_w_dry,\
     101from anuga.config import g, epsilon, beta_h, beta_w, beta_w_dry,\
    102102     beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
    103103from anuga.config import alpha_balance
    … …  
    563563    """
    564564
    565     from anuga.config import g, epsilon
    566565    from math import sqrt
    567566
    … …  
    657656    """
    658657
    659     from anuga.config import g, epsilon
    660658    from math import sqrt
    661659    from Numeric import array
    … …  
    835833    N = len(domain) # number_of_triangles
    836834
    837     #Shortcuts
     835    # Shortcuts
    838836    Stage = domain.quantities['stage']
    839837    Xmom = domain.quantities['xmomentum']
    840838    Ymom = domain.quantities['ymomentum']
    841839    Elevation = domain.quantities['elevation']
     840
    842841    from shallow_water_ext import extrapolate_second_order_sw
    843842    extrapolate_second_order_sw(domain,
    … …  
    853852                                Xmom.vertex_values,
    854853                                Ymom.vertex_values,
    855                                 Elevation.vertex_values)
     854                                Elevation.vertex_values,
     855                                int(domain.optimise_dry_cells))
    856856
    857857def compute_fluxes_c(domain):

  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4690 r4710  
    4848}
    4949
    50 int find_qmin_and_qmax(double dq0, double dq1, double dq2, double *qmin, double *qmax){
    51   //Considering the centroid of an FV triangle and the vertices of its auxiliary triangle, find
    52   //qmin=min(q)-qc and qmax=max(q)-qc, where min(q) and max(q) are respectively min and max over the
    53   //four values (at the centroid of the FV triangle and the auxiliary triangle vertices),
    54   //and qc is the centroid
    55   //dq0=q(vertex0)-q(centroid of FV triangle)
    56   //dq1=q(vertex1)-q(vertex0)
    57   //dq2=q(vertex2)-q(vertex0)
     50int find_qmin_and_qmax(double dq0, double dq1, double dq2,
     51                       double *qmin, double *qmax){
     52  // Considering the centroid of an FV triangle and the vertices of its
     53  // auxiliary triangle, find
     54  // qmin=min(q)-qc and qmax=max(q)-qc,
     55  // where min(q) and max(q) are respectively min and max over the
     56  // four values (at the centroid of the FV triangle and the auxiliary
     57  // triangle vertices),
     58  // and qc is the centroid
     59  // dq0=q(vertex0)-q(centroid of FV triangle)
     60  // dq1=q(vertex1)-q(vertex0)
     61  // dq2=q(vertex2)-q(vertex0)
    5862  if (dq0>=0.0){
    5963    if (dq1>=dq2){
    … …  
    6367        *qmax=dq0;
    6468      if ((*qmin=dq0+dq2)<0)
    65         ;//qmin is already set to correct value
     69        ;// qmin is already set to correct value
    6670      else
    6771        *qmin=0.0;
    6872    }
    69     else{//dq1<dq2
     73    else{// dq1<dq2
    7074      if (dq2>0)
    7175        *qmax=dq0+dq2;
    … …  
    7377        *qmax=dq0;
    7478      if ((*qmin=dq0+dq1)<0)
    75         ;//qmin is the correct value
     79        ;// qmin is the correct value
    7680      else
    7781        *qmin=0.0;
    … …  
    8589        *qmin=dq0;
    8690      if ((*qmax=dq0+dq2)>0.0)
    87         ;//qmax is already set to the correct value
     91        ;// qmax is already set to the correct value
    8892      else
    8993        *qmax=0.0;
    9094    }
    91     else{//dq1>dq2
     95    else{// dq1>dq2
    9296      if (dq2<0.0)
    9397        *qmin=dq0+dq2;
    … …  
    9599        *qmin=dq0;
    96100      if ((*qmax=dq0+dq1)>0.0)
    97         ;//qmax is already set to the correct value
     101        ;// qmax is already set to the correct value
    98102      else
    99103        *qmax=0.0;
    … …  
    104108
    105109int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
    106   //given provisional jumps dqv from the FV triangle centroid to its vertices and
    107   //jumps qmin (qmax) between the centroid of the FV triangle and the
    108   //minimum (maximum) of the values at the centroid of the FV triangle and the auxiliary triangle vertices,
    109   //calculate a multiplicative factor phi by which the provisional vertex jumps are to be limited
     110  // Given provisional jumps dqv from the FV triangle centroid to its
     111  // vertices and jumps qmin (qmax) between the centroid of the FV
     112  // triangle and the minimum (maximum) of the values at the centroid of
     113  // the FV triangle and the auxiliary triangle vertices,
     114  // calculate a multiplicative factor phi by which the provisional
     115  // vertex jumps are to be limited
    110116  int i;
    111117  double r=1000.0, r0=1.0, phi=1.0;
    112   static double TINY = 1.0e-100;//to avoid machine accuracy problems.
    113   //Any provisional jump with magnitude < TINY does not contribute to the limiting process.
    114   for (i=0;i<3;i++){
     118  static double TINY = 1.0e-100; // to avoid machine accuracy problems.
     119  // FIXME: Perhaps use the epsilon used elsewhere.
     120 
     121  // Any provisional jump with magnitude < TINY does not contribute to
     122  // the limiting process.
     123    for (i=0;i<3;i++){
    115124    if (dqv[i]<-TINY)
    116125      r0=qmin/dqv[i];
    … …  
    118127      r0=qmax/dqv[i];
    119128    r=min(r0,r);
    120     //
    121   }
     129  }
     130 
    122131  phi=min(r*beta_w,1.0);
    123132  for (i=0;i<3;i++)
    … …  
    866875  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
    867876  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
    868   double dqv[3], qmin, qmax, hmin;
     877  double dqv[3], qmin, qmax, hmin, hmax;
    869878  double hc, h0, h1, h2;
    870   double epsilon=1.0e-12; // FIXME Pass in
     879  double epsilon=1.0e-12;
     880  int optimise_dry_cells=1; // Optimisation flag   
    871881  double beta_w, beta_w_dry, beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, beta_tmp;
    872882  double minimum_allowed_height;
    873   //provisional jumps from centroids to v'tices and safety factor re limiting
    874   //by which these jumps are limited
     883  // Provisional jumps from centroids to v'tices and safety factor re limiting
     884  // by which these jumps are limited
    875885  // Convert Python arguments to C
    876   if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOO",
     886  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
    877887                        &domain,
    878888                        &surrogate_neighbours,
    … …  
    887897                        &xmom_vertex_values,
    888898                        &ymom_vertex_values,
    889                         &elevation_vertex_values)) {
    890     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    891     return NULL;
    892   }
    893 
    894   //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
     899                        &elevation_vertex_values,
     900                        &optimise_dry_cells)) {                 
     901                       
     902    PyErr_SetString(PyExc_RuntimeError, "Input arguments to extrapolate_second_order_sw failed");
     903    return NULL;
     904  }
     905
     906  // FIXME (Ole): Investigate if it is quicker to obtain all input arguments using GetAttrString rather than ParseTuple.
     907  // It certainly looked as if passing domain.epsilon is slowed things down
     908 
     909  // Get the safety factor beta_w, set in the config.py file. This is used in the limiting process
    895910  Tmp = PyObject_GetAttrString(domain, "beta_w");
    896911  if (!Tmp) {
    … …  
    943958  Tmp = PyObject_GetAttrString(domain, "minimum_allowed_height");
    944959  if (!Tmp) {
    945     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_heigt");
     960    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object minimum_allowed_height");
    946961    return NULL;
    947962  } 
    948963  minimum_allowed_height = PyFloat_AsDouble(Tmp);
    949964  Py_DECREF(Tmp); 
    950  
     965
     966  Tmp = PyObject_GetAttrString(domain, "epsilon");
     967  if (!Tmp) {
     968    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: extrapolate_second_order_sw could not obtain object epsilon");
     969    return NULL;
     970  } 
     971  epsilon = PyFloat_AsDouble(Tmp);
     972  Py_DECREF(Tmp); 
     973   
    951974  number_of_elements = stage_centroid_values -> dimensions[0];
    952975  for (k=0; k<number_of_elements; k++) {
    … …  
    10341057      h1 = ((double *)stage_centroid_values->data)[k1] - ((double *)elevation_centroid_values->data)[k1];
    10351058      h2 = ((double *)stage_centroid_values->data)[k2] - ((double *)elevation_centroid_values->data)[k2];
    1036       hmin = min(hc,min(h0,min(h1,h2)));
    1037      
    1038      
    1039       // Dry cell optimisation experiment
    1040       //printf("hmin = %e\n", hmin);     
    1041       //if (hmin < epsilon) {
    1042         //printf("Dry\n");
    1043         //continue;
    1044       //}
     1059      hmin = min(hc,min(h0,min(h1,h2)));  // FIXME Don't need to include hc
     1060     
     1061     
     1062      if (optimise_dry_cells) {     
     1063        // Check if linear reconstruction is necessary for triangle k
     1064        // This check will exclude dry cells.
     1065
     1066        hmax = max(h0,max(h1,h2));     
     1067        if (hmax < epsilon) {
     1068          continue;
     1069        }
     1070      }
    10451071
    10461072           

  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r4704 r4710  
    28592859
    28602860        domain.distribute_to_vertices_and_edges()
     2861
    28612862        assert allclose(domain.quantities['stage'].vertex_values[:12,0],
    28622863                        [0.0001714, 0.02656103, 0.00024152,
    … …  
    49834984if __name__ == "__main__":
    49844985
    4985     #suite = unittest.makeSuite(Test_Shallow_Water,'test')   
    4986     suite = unittest.makeSuite(Test_Shallow_Water,'test_extrema')   
     4986    suite = unittest.makeSuite(Test_Shallow_Water,'test')   
     4987    #suite = unittest.makeSuite(Test_Shallow_Water,'test_extrema')   
    49874988    #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
    49884989    #suite = unittest.makeSuite(Test_Shallow_Water,'test_get_maximum_inundation_from_sww')
Note: See TracChangeset for help on using the changeset viewer.