Changeset 1641

Timestamp:

Jul 26, 2005, 5:16:04 PM (20 years ago)

Author:

ole

Message:

Undid changes committed in r1458 and r1459:

svn merge -r 1459:1457 ​https://datamining/...../pyvolution

and thus brought back Matt's optimisations

modified shallow_water.py
modified shallow_water_ext.c

Location:

inundation/ga/storm_surge/pyvolution

Files:

• shallow_water.py (modified) (13 diffs)
• shallow_water_ext.c (modified) (16 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -56 +56 @@
 
 from domain import *
-from region import *
+from region import *#
+
 Generic_domain = Domain #Rename
 
@@ -81 +82 @@
         #Realtime visualisation
         self.visualiser = None
-        self.visualise = False
+        self.visualise  = False
         self.visualise_color_stage = False
         self.visualise_stage_range = 1.0
@@ -107 +108 @@
         #self.eta = self.quantities['friction']
 
+    def initialise_visualiser(self,scale_z=1.0,rect=None):
+        #Realtime visualisation
+        if self.visualiser is None:
+            from realtime_visualisation_new import Visualiser
+            self.visualiser = Visualiser(self,scale_z,rect)
+        self.visualise = True
+
     def check_integrity(self):
         Generic_domain.check_integrity(self)
@@ -119 +127 @@
         assert self.conserved_quantities[2] == 'ymomentum', msg
 
+    def extrapolate_second_order_sw(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order_sw(self)
 
     def compute_fluxes(self):
@@ -250 +262 @@
         self.writer.store_timestep(name)
 
+####################MH 090605 new extrapolation function belonging to domain class
+def extrapolate_second_order_sw(domain):
+    """extrapolate conserved quntities to the vertices of the triangles
+    Python version to be written after the C version
+    """
+    msg = 'Method extrapolate_second_order_sw should be implemented in C'
+    raise msg
+####################MH 090605 ###########################################
 
     def initialise_visualiser(self,scale_z=1.0,rect=None):
@@ -436 +456 @@
 
     flux = zeros(3, Float) #Work array for summing up fluxes
+
 
     #Loop
@@ -482 +503 @@
     domain.timestep = timestep
 
+#MH090605 The following method belongs to the shallow_water domain class
+#see comments in the corresponding method in shallow_water_ext.c
+def extrapolate_second_order_sw_c(domain):
+    """Wrapper calling C version of extrapolate_second_order_sw
+    """
+    import sys
+    from Numeric import zeros, Float
+
+    N = domain.number_of_elements
+
+    #Shortcuts
+    Stage = domain.quantities['stage']
+    Xmom = domain.quantities['xmomentum']
+    Ymom = domain.quantities['ymomentum']
+    from shallow_water_ext import extrapolate_second_order_sw
+    extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
+                                domain.number_of_boundaries,
+                                domain.centroid_coordinates,
+                                Stage.centroid_values,
+                                Xmom.centroid_values,
+                                Ymom.centroid_values,
+                                domain.vertex_coordinates,
+                                Stage.vertex_values,
+                                Xmom.vertex_values,
+                                Ymom.vertex_values)
 
 def compute_fluxes_c(domain):
@@ -516 +562 @@
                                      Stage.explicit_update,
                                      Xmom.explicit_update,
-                                     Ymom.explicit_update)
+                                     Ymom.explicit_update,
+                                     domain.already_computed_flux)
 
 
@@ -545 +592 @@
     """
 
+    from config import optimised_gradient_limiter
+
     #Remove very thin layers of water
     protect_against_infinitesimal_and_negative_heights(domain)
 
     #Extrapolate all conserved quantities
-    for name in domain.conserved_quantities:
-        Q = domain.quantities[name]
-        if domain.order == 1:
-            Q.extrapolate_first_order()
+    if optimised_gradient_limiter:
+        #MH090605 if second order,
+        #perform the extrapolation and limiting on
+        #all of the conserved quantitie
+        
+        if (domain.order == 1):
+            for name in domain.conserved_quantities:
+                Q = domain.quantities[name]
+                Q.extrapolate_first_order()
         elif domain.order == 2:
-            Q.extrapolate_second_order()
-            Q.limit()
+            domain.extrapolate_second_order_sw()
         else:
             raise 'Unknown order'
+    else:    
+        #old code:
+        for name in domain.conserved_quantities:
+            Q = domain.quantities[name]
+            if domain.order == 1:
+                Q.extrapolate_first_order()
+            elif domain.order == 2:
+                Q.extrapolate_second_order()
+                Q.limit()
+            else:
+                raise 'Unknown order'
+            
 
     #Take bed elevation into account when water heights are small
@@ -738 +803 @@
 
     #Call C-extension
-    from shallow_water_ext import h_limiter
+    from shallow_water_ext import h_limiter_sw as h_limiter
     hvbar = h_limiter(domain, hc, hv)
 
@@ -909 +974 @@
             raise msg
 
-        #Handy shorthands
-        self.stage = domain.quantities['stage'].edge_values
-        self.xmom = domain.quantities['xmomentum'].edge_values
-        self.ymom = domain.quantities['ymomentum'].edge_values
-        self.normals = domain.normals
-
-        from Numeric import zeros, Float
-        self.conserved_quantities = zeros(3, Float)
+        #Handy shorthands
+        self.stage   = domain.quantities['stage'].edge_values
+        self.xmom    = domain.quantities['xmomentum'].edge_values
+        self.ymom    = domain.quantities['ymomentum'].edge_values
+        self.normals = domain.normals
+
+        from Numeric import zeros, Float
+        self.conserved_quantities = zeros(3, Float)
 
     def __repr__(self):
@@ -927 +992 @@
         """
 
-        q = self.conserved_quantities
-        q[0] = self.stage[vol_id, edge_id]
-        q[1] = self.xmom[vol_id, edge_id]
-        q[2] = self.ymom[vol_id, edge_id]
-
-        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
+        q = self.conserved_quantities
+        q[0] = self.stage[vol_id, edge_id]
+        q[1] = self.xmom[vol_id, edge_id]
+        q[2] = self.ymom[vol_id, edge_id]
+
+        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
 
 
@@ -1647 +1712 @@
     from shallow_water_ext import rotate, assign_windfield_values
     compute_fluxes = compute_fluxes_c
+    extrapolate_second_order_sw=extrapolate_second_order_sw_c
     gravity = gravity_c
     manning_friction = manning_friction_c

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -47 +47 @@
 }
 
-
+int find_qmin_and_qmax(double dq0, double dq1, double dq2, double *qmin, double *qmax){
+  //Considering the centroid of an FV triangle and the vertices of its auxiliary triangle, find
+  //qmin=min(q)-qc and qmax=max(q)-qc, where min(q) and max(q) are respectively min and max over the
+  //four values (at the centroid of the FV triangle and the auxiliary triangle vertices),
+  //and qc is the centroid
+  //dq0=q(vertex0)-q(centroid of FV triangle)
+  //dq1=q(vertex1)-q(vertex0)
+  //dq2=q(vertex2)-q(vertex0)
+  if (dq0>=0.0){
+    if (dq1>=dq2){
+      if (dq1>=0.0)
+        *qmax=dq0+dq1;
+      else
+        *qmax=dq0;
+      if ((*qmin=dq0+dq2)<0)
+        ;//qmin is already set to correct value
+      else
+        *qmin=0.0;
+    }
+    else{//dq1<dq2
+      if (dq2>0)
+        *qmax=dq0+dq2;
+      else
+        *qmax=dq0;
+      if ((*qmin=dq0+dq1)<0)
+        ;//qmin is the correct value
+      else
+        *qmin=0.0;
+    }
+  }
+  else{//dq0<0
+    if (dq1<=dq2){
+      if (dq1<0.0)
+        *qmin=dq0+dq1;
+      else
+        *qmin=dq0;
+      if ((*qmax=dq0+dq2)>0.0)
+        ;//qmax is already set to the correct value
+      else
+        *qmax=0.0;
+    }
+    else{//dq1>dq2
+      if (dq2<0.0)
+        *qmin=dq0+dq2;
+      else
+        *qmin=dq0;
+      if ((*qmax=dq0+dq1)>0.0)
+        ;//qmax is already set to the correct value
+      else
+        *qmax=0.0;
+    }
+  }
+  return 0;
+}
+
+int limit_gradient(double *dqv, double qmin, double qmax, double beta_w){
+  //given provisional jumps dqv from the FV triangle centroid to its vertices and
+  //jumps qmin (qmax) between the centroid of the FV triangle and the
+  //minimum (maximum) of the values at the centroid of the FV triangle and the auxiliary triangle vertices,
+  //calculate a multiplicative factor phi by which the provisional vertex jumps are to be limited
+  int i;
+  double r=1000.0, r0=1.0, phi=1.0;
+  static double TINY = 1.0e-100;//to avoid machine accuracy problems.
+  //Any provisional jump with magnitude < TINY does not contribute to the limiting process.
+  for (i=0;i<3;i++){
+    if (dqv[i]<-TINY)
+      r0=qmin/dqv[i];
+    if (dqv[i]>TINY)
+      r0=qmax/dqv[i];
+    r=min(r0,r);
+    //
+  }
+  phi=min(r*beta_w,1.0);
+  for (i=0;i<3;i++)
+    dqv[i]=dqv[i]*phi;
+  return 0;
+}
 
 // Computational function for flux computation (using stage w=z+h)
@@ -170 +246 @@
       if (h >= eps) {
         S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
-        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
+        //S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
+        S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
         //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
 
@@ -423 +500 @@
 }
 
+PyObject *extrapolate_second_order_sw(PyObject *self, PyObject *args) {
+  /*Compute the vertex values based on a linear reconstruction on each triangle
+    These values are calculated as follows:
+    1) For each triangle not adjacent to a boundary, we consider the auxiliary triangle
+    formed by the centroids of its three neighbours.
+    2) For each conserved quantity, we integrate around the auxiliary triangle's boundary the product
+    of the quantity and the outward normal vector. Dividing by the triangle area gives (a,b), the average
+    of the vector (q_x,q_y) on the auxiliary triangle. We suppose that the linear reconstruction on the
+    original triangle has gradient (a,b).
+    3) Provisional vertex junmps dqv[0,1,2] are computed and these are then limited by calling the functions
+    find_qmin_and_qmax and limit_gradient
+
+    Python call:
+    extrapolate_second_order_sw(domain.surrogate_neighbours,
+                                domain.number_of_boundaries
+                                domain.centroid_coordinates,
+                                Stage.centroid_values
+                                Xmom.centroid_values
+                                Ymom.centroid_values
+                                domain.vertex_coordinates,
+                                Stage.vertex_values,
+                                Xmom.vertex_values,
+                                Ymom.vertex_values)
+
+    Post conditions:
+            The vertices of each triangle have values from a limited linear reconstruction
+            based on centroid values
+
+  */
+  PyArrayObject *surrogate_neighbours,
+    *number_of_boundaries,
+    *centroid_coordinates,
+    *stage_centroid_values,
+    *xmom_centroid_values,
+    *ymom_centroid_values,
+    *vertex_coordinates,
+    *stage_vertex_values,
+    *xmom_vertex_values,
+    *ymom_vertex_values;
+  PyObject *domain, *Tmp;
+  //Local variables
+  double a, b;//gradient vector, not stored but used to calculate vertex values from centroids
+  int number_of_elements,k,k0,k1,k2,k3,k6,coord_index,i;
+  double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2;//vertices of the auxiliary triangle
+  double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
+  double dqv[3], qmin, qmax, beta_w;//provisional jumps from centroids to v'tices and safety factor re limiting
+  //by which these jumps are limited
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO",
+                        &domain,
+                        &surrogate_neighbours,
+                        &number_of_boundaries,
+                        &centroid_coordinates,
+                        &stage_centroid_values,
+                        &xmom_centroid_values,
+                        &ymom_centroid_values,
+                        &vertex_coordinates,
+                        &stage_vertex_values,
+                        &xmom_vertex_values,
+                        &ymom_vertex_values)) {
+    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+    return NULL;
+  }
+
+  //get the safety factor beta_w, set in the config.py file. This is used in the limiting process
+  Tmp = PyObject_GetAttrString(domain, "beta_w");
+  if (!Tmp)
+    return NULL;
+  beta_w = PyFloat_AsDouble(Tmp);
+  Py_DECREF(Tmp);
+  number_of_elements = stage_centroid_values -> dimensions[0];
+  for (k=0; k<number_of_elements; k++) {
+    k3=k*3;
+    k6=k*6;
+
+    if (((long *) number_of_boundaries->data)[k]==3){/*no neighbours, set gradient on the triangle to zero*/
+      ((double *) stage_vertex_values->data)[k3]=((double *)stage_centroid_values->data)[k];
+      ((double *) stage_vertex_values->data)[k3+1]=((double *)stage_centroid_values->data)[k];
+      ((double *) stage_vertex_values->data)[k3+2]=((double *)stage_centroid_values->data)[k];
+      ((double *) xmom_vertex_values->data)[k3]=((double *)xmom_centroid_values->data)[k];
+      ((double *) xmom_vertex_values->data)[k3+1]=((double *)xmom_centroid_values->data)[k];
+      ((double *) xmom_vertex_values->data)[k3+2]=((double *)xmom_centroid_values->data)[k];
+      ((double *) ymom_vertex_values->data)[k3]=((double *)ymom_centroid_values->data)[k];
+      ((double *) ymom_vertex_values->data)[k3+1]=((double *)ymom_centroid_values->data)[k];
+      ((double *) ymom_vertex_values->data)[k3+2]=((double *)ymom_centroid_values->data)[k];
+      continue;
+    }
+    else {//we will need centroid coordinates and vertex coordinates of the triangle
+      //get the vertex coordinates of the FV triangle
+      xv0=((double *)vertex_coordinates->data)[k6]; yv0=((double *)vertex_coordinates->data)[k6+1];
+      xv1=((double *)vertex_coordinates->data)[k6+2]; yv1=((double *)vertex_coordinates->data)[k6+3];
+      xv2=((double *)vertex_coordinates->data)[k6+4]; yv2=((double *)vertex_coordinates->data)[k6+5];
+      //get the centroid coordinates of the FV triangle
+      coord_index=2*k;
+      x=((double *)centroid_coordinates->data)[coord_index];
+      y=((double *)centroid_coordinates->data)[coord_index+1];
+      //store x- and y- differentials for the vertices of the FV triangle relative to the centroid
+      dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
+      dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
+    }
+    if (((long *)number_of_boundaries->data)[k]<=1){
+      //if no boundaries, auxiliary triangle is formed from the centroids of the three neighbours
+      //if one boundary, auxiliary triangle is formed from this centroid and its two neighbours
+      k0=((long *)surrogate_neighbours->data)[k3];
+      k1=((long *)surrogate_neighbours->data)[k3+1];
+      k2=((long *)surrogate_neighbours->data)[k3+2];
+      //get the auxiliary triangle's vertex coordinates (really the centroids of neighbouring triangles)
+      coord_index=2*k0;
+      x0=((double *)centroid_coordinates->data)[coord_index];
+      y0=((double *)centroid_coordinates->data)[coord_index+1];
+      coord_index=2*k1;
+      x1=((double *)centroid_coordinates->data)[coord_index];
+      y1=((double *)centroid_coordinates->data)[coord_index+1];
+      coord_index=2*k2;
+      x2=((double *)centroid_coordinates->data)[coord_index];
+      y2=((double *)centroid_coordinates->data)[coord_index+1];
+      //store x- and y- differentials for the vertices of the auxiliary triangle
+      dx1=x1-x0; dx2=x2-x0;
+      dy1=y1-y0; dy2=y2-y0;
+      //calculate 2*area of the auxiliary triangle
+      area2 = dy2*dx1 - dy1*dx2;//the triangle is guaranteed to be counter-clockwise
+      //If the mesh is 'weird' near the boundary, the trianlge might be flat or clockwise:
+      if (area2<=0)
+        return NULL;
+
+      //### stage ###
+      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
+      dq0=((double *)stage_centroid_values->data)[k0]-((double *)stage_centroid_values->data)[k];
+      //calculate differentials between the vertices of the auxiliary triangle
+      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k0];
+      dq2=((double *)stage_centroid_values->data)[k2]-((double *)stage_centroid_values->data)[k0];
+      //calculate the gradient of stage on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a /= area2;
+      b = dx1*dq2 - dx2*dq1;
+      b /= area2;
+      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
+      dqv[0]=a*dxv0+b*dyv0;
+      dqv[1]=a*dxv1+b*dyv1;
+      dqv[2]=a*dxv2+b*dyv2;
+      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
+      //and compute jumps from the centroid to the min and max
+      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
+      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
+
+      //### xmom ###
+      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
+      dq0=((double *)xmom_centroid_values->data)[k0]-((double *)xmom_centroid_values->data)[k];
+      //calculate differentials between the vertices of the auxiliary triangle
+      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k0];
+      dq2=((double *)xmom_centroid_values->data)[k2]-((double *)xmom_centroid_values->data)[k0];
+      //calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a /= area2;
+      b = dx1*dq2 - dx2*dq1;
+      b /= area2;
+      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
+      dqv[0]=a*dxv0+b*dyv0;
+      dqv[1]=a*dxv1+b*dyv1;
+      dqv[2]=a*dxv2+b*dyv2;
+      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
+      //and compute jumps from the centroid to the min and max
+      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
+      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
+
+      //### ymom ###
+      //calculate the difference between vertex 0 of the auxiliary triangle and the FV triangle centroid
+      dq0=((double *)ymom_centroid_values->data)[k0]-((double *)ymom_centroid_values->data)[k];
+      //calculate differentials between the vertices of the auxiliary triangle
+      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k0];
+      dq2=((double *)ymom_centroid_values->data)[k2]-((double *)ymom_centroid_values->data)[k0];
+      //calculate the gradient of xmom on the auxiliary triangle
+      a = dy2*dq1 - dy1*dq2;
+      a /= area2;
+      b = dx1*dq2 - dx2*dq1;
+      b /= area2;
+      //calculate provisional jumps in stage from the centroid of the FV tri to its vertices, to be limited
+      dqv[0]=a*dxv0+b*dyv0;
+      dqv[1]=a*dxv1+b*dyv1;
+      dqv[2]=a*dxv2+b*dyv2;
+      //now we want to find min and max of the centroid and the vertices of the auxiliary triangle
+      //and compute jumps from the centroid to the min and max
+      find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
+      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
+    }//if (number_of_boundaries[k]<=1)
+    else{//number_of_boundaries==2
+      //one internal neighbour and gradient is in direction of the neighbour's centroid
+      //find the only internal neighbour
+      for (k2=k3;k2<k3+3;k2++){//k2 just indexes the edges of triangle k
+        if (((long *)surrogate_neighbours->data)[k2]!=k)//find internal neighbour of triabngle k
+          break;
+      }
+      if ((k2==k3+3))//if we didn't find an internal neighbour
+        return NULL;//error
+      k1=((long *)surrogate_neighbours->data)[k2];
+      //the coordinates of the triangle are already (x,y). Get centroid of the neighbour (x1,y1)
+      coord_index=2*k1;
+      x1=((double *)centroid_coordinates->data)[coord_index];
+      y1=((double *)centroid_coordinates->data)[coord_index+1];
+      //compute x- and y- distances between the centroid of the FV triangle and that of its neighbour
+      dx1=x1-x; dy1=y1-y;
+      //set area2 as the square of the distance
+      area2=dx1*dx1+dy1*dy1;
+      //set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2) and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
+      //respectively correspond to the x- and y- gradients of the conserved quantities
+      dx2=1.0/area2;
+      dy2=dx2*dy1;
+      dx2*=dx1;
+
+      //## stage ###
+      //compute differentials
+      dq1=((double *)stage_centroid_values->data)[k1]-((double *)stage_centroid_values->data)[k];
+      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
+      a=dq1*dx2;
+      b=dq1*dy2;
+      //calculate provisional vertex jumps, to be limited
+      dqv[0]=a*dxv0+b*dyv0;
+      dqv[1]=a*dxv1+b*dyv1;
+      dqv[2]=a*dxv2+b*dyv2;
+      //now limit the jumps
+      if (dq1>=0.0){
+        qmin=0.0;
+        qmax=dq1;
+      }
+      else{
+        qmin=dq1;
+        qmax=0.0;
+      }
+      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+        ((double *)stage_vertex_values->data)[k3+i]=((double *)stage_centroid_values->data)[k]+dqv[i];
+
+      //## xmom ###
+      //compute differentials
+      dq1=((double *)xmom_centroid_values->data)[k1]-((double *)xmom_centroid_values->data)[k];
+      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
+      a=dq1*dx2;
+      b=dq1*dy2;
+      //calculate provisional vertex jumps, to be limited
+      dqv[0]=a*dxv0+b*dyv0;
+      dqv[1]=a*dxv1+b*dyv1;
+      dqv[2]=a*dxv2+b*dyv2;
+      //now limit the jumps
+      if (dq1>=0.0){
+        qmin=0.0;
+        qmax=dq1;
+      }
+      else{
+        qmin=dq1;
+        qmax=0.0;
+      }
+      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+        ((double *)xmom_vertex_values->data)[k3+i]=((double *)xmom_centroid_values->data)[k]+dqv[i];
+
+      //## ymom ###
+      //compute differentials
+      dq1=((double *)ymom_centroid_values->data)[k1]-((double *)ymom_centroid_values->data)[k];
+      //calculate the gradient between the centroid of the FV triangle and that of its neighbour
+      a=dq1*dx2;
+      b=dq1*dy2;
+      //calculate provisional vertex jumps, to be limited
+      dqv[0]=a*dxv0+b*dyv0;
+      dqv[1]=a*dxv1+b*dyv1;
+      dqv[2]=a*dxv2+b*dyv2;
+      //now limit the jumps
+      if (dq1>=0.0){
+        qmin=0.0;
+        qmax=dq1;
+      }
+      else{
+        qmin=dq1;
+        qmax=0.0;
+      }
+      limit_gradient(dqv,qmin,qmax,beta_w);//the gradient will be limited
+      for (i=0;i<3;i++)
+        ((double *)ymom_vertex_values->data)[k3+i]=((double *)ymom_centroid_values->data)[k]+dqv[i];
+    }//else [number_of_boudaries==2]
+  }//for k=0 to number_of_elements-1
+  return Py_BuildValue("");
+}//extrapolate_second-order_sw
+
 PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
   //
@@ -480 +845 @@
 }
 
-
 PyObject *compute_fluxes(PyObject *self, PyObject *args) {
   /*Compute all fluxes and the timestep suitable for all volumes
@@ -499 +863 @@
     domain.timestep = compute_fluxes(timestep,
                                      domain.epsilon,
-                                     domain.g,
+                                     domain.g,
                                      domain.neighbours,
                                      domain.neighbour_edges,
@@ -515 +879 @@
                                      Stage.explicit_update,
                                      Xmom.explicit_update,
-                                     Ymom.explicit_update)
+                                     Ymom.explicit_update,
+                                     already_computed_flux)
 
 
@@ -537 +902 @@
     *stage_explicit_update,
     *xmom_explicit_update,
-    *ymom_explicit_update;
+    *ymom_explicit_update,
+    *already_computed_flux;//tracks whether the flux across an edge has already been computed
 
 
@@ -543 +909 @@
   double timestep, max_speed, epsilon, g;
   double normal[2], ql[3], qr[3], zl, zr;
-  double flux[3], edgeflux[3]; //Work arrays for summing up fluxes
-
-  int number_of_elements, k, i, j, m, n;
-  int ki, nm, ki2; //Index shorthands
+  double edgeflux[3]; //Work arrays for summing up fluxes
+
+  int number_of_elements, k, i, m, n;
+  int ki, nm=0, ki2; //Index shorthands
+  static long call=1;
 
 
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOO",
+  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOO",
                         &timestep,
                         &epsilon,
@@ -567 +934 @@
                         &stage_explicit_update,
                         &xmom_explicit_update,
-                        &ymom_explicit_update)) {
+                        &ymom_explicit_update,
+                        &already_computed_flux)) {
     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     return NULL;
   }
-
   number_of_elements = stage_edge_values -> dimensions[0];
-
-
+  call++;//a static local variable to which already_computed_flux is compared
+  //set explicit_update to zero for all conserved_quantities.
+  //This assumes compute_fluxes called before forcing terms
   for (k=0; k<number_of_elements; k++) {
-
-    //Reset work array
-    for (j=0; j<3; j++) flux[j] = 0.0;
-
-    //Loop through neighbours and compute edge flux for each
+    ((double *) stage_explicit_update -> data)[k]=0.0;
+    ((double *) xmom_explicit_update -> data)[k]=0.0;
+    ((double *) ymom_explicit_update -> data)[k]=0.0;
+  }
+  //Loop through neighbours and compute edge flux for each
+  for (k=0; k<number_of_elements; k++) {
     for (i=0; i<3; i++) {
       ki = k*3+i;
+      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
+        continue;
       ql[0] = ((double *) stage_edge_values -> data)[ki];
       ql[1] = ((double *) xmom_edge_values -> data)[ki];
@@ -598 +969 @@
       } else {
         m = ((long *) neighbour_edges -> data)[ki];
-
         nm = n*3+m;
         qr[0] = ((double *) stage_edge_values -> data)[nm];
@@ -605 +975 @@
         zr =    ((double *) bed_edge_values -> data)[nm];
       }
-
-      //printf("%d %d [%d] (%d, %d): %.2f %.2f %.2f\n", k, i, ki, n, m,
-      //             ql[0], ql[1], ql[2]);
-
-
       // Outward pointing normal vector
       // normal = domain.normals[k, 2*i:2*i+2]
@@ -615 +980 @@
       normal[0] = ((double *) normals -> data)[ki2];
       normal[1] = ((double *) normals -> data)[ki2+1];
-
       //Edge flux computation
       flux_function(ql, qr, zl, zr,
@@ -621 +985 @@
                     epsilon, g,
                     edgeflux, &max_speed);
-
-
-      //flux -= edgeflux * edgelengths[k,i]
-      for (j=0; j<3; j++) {
-        flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
+      //update triangle k
+      ((long *) already_computed_flux->data)[ki]=call;
+      ((double *) stage_explicit_update -> data)[k] -= edgeflux[0]*((double *) edgelengths -> data)[ki];
+      ((double *) xmom_explicit_update -> data)[k] -= edgeflux[1]*((double *) edgelengths -> data)[ki];
+      ((double *) ymom_explicit_update -> data)[k] -= edgeflux[2]*((double *) edgelengths -> data)[ki];
+      //update the neighbour n
+      if (n>=0){
+        ((long *) already_computed_flux->data)[nm]=call;
+        ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
+        ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
+        ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
       }
-
-      //Update timestep
-      //timestep = min(timestep, domain.radii[k]/max_speed)
-      //FIXME: SR Add parameter for CFL condition
-      if (max_speed > epsilon) {
-        timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
-      }
+      ///for (j=0; j<3; j++) {
+        ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
+        ///}
+        //Update timestep
+        //timestep = min(timestep, domain.radii[k]/max_speed)
+        //FIXME: SR Add parameter for CFL condition
+        if (max_speed > epsilon) {
+          timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
+          //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
+          if (n>=0)
+            timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
+        }
     } // end for i
-
     //Normalise by area and store for when all conserved
     //quantities get updated
-    // flux /= areas[k]
-    for (j=0; j<3; j++) {
-      flux[j] /= ((double *) areas -> data)[k];
-    }
-
-    ((double *) stage_explicit_update -> data)[k] = flux[0];
-    ((double *) xmom_explicit_update -> data)[k] = flux[1];
-    ((double *) ymom_explicit_update -> data)[k] = flux[2];
-
+    ((double *) stage_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
+    ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
+    ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
   } //end for k
-
   return Py_BuildValue("d", timestep);
 }
-
-
 
 PyObject *protect(PyObject *self, PyObject *args) {
@@ -812 +1177 @@
 }
 
-
-
+PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
+  //a faster version of h_limiter above
+  PyObject *domain, *Tmp;
+  PyArrayObject
+    *hv, *hc, //Depth at vertices and centroids
+    *hvbar,   //Limited depth at vertices (return values)
+    *neighbours;
+
+  int k, i, N, k3,k0,k1,k2;
+  int dimensions[2];
+  double beta_h; //Safety factor (see config.py)
+  double hmin, hmax, dh[3];
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
+    return NULL;
+  neighbours = get_consecutive_array(domain, "neighbours");
+
+  //Get safety factor beta_h
+  Tmp = PyObject_GetAttrString(domain, "beta_h");
+  if (!Tmp)
+    return NULL;
+  beta_h = PyFloat_AsDouble(Tmp);
+
+  Py_DECREF(Tmp);
+
+  N = hc -> dimensions[0];
+
+  //Create hvbar
+  dimensions[0] = N;
+  dimensions[1] = 3;
+  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
+  for (k=0;k<N;k++){
+    k3=k*3;
+    //get the ids of the neighbours
+    k0 = ((long*) neighbours -> data)[k3];
+    k1 = ((long*) neighbours -> data)[k3+1];
+    k2 = ((long*) neighbours -> data)[k3+2];
+    //set hvbar provisionally
+    for (i=0;i<3;i++){
+      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
+      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
+    }
+    hmin=((double*) hc -> data)[k];
+    hmax=hmin;
+    if (k0>=0){
+      hmin=min(hmin,((double*) hc -> data)[k0]);
+      hmax=max(hmax,((double*) hc -> data)[k0]);
+    }
+    if (k1>=0){
+      hmin=min(hmin,((double*) hc -> data)[k1]);
+      hmax=max(hmax,((double*) hc -> data)[k1]);
+    }
+    if (k2>=0){
+      hmin=min(hmin,((double*) hc -> data)[k2]);
+      hmax=max(hmax,((double*) hc -> data)[k2]);
+    }
+    hmin-=((double*) hc -> data)[k];
+    hmax-=((double*) hc -> data)[k];
+    limit_gradient(dh,hmin,hmax,beta_h);
+    for (i=0;i<3;i++)
+      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
+  }
+  return PyArray_Return(hvbar);
+}
 
 PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
@@ -864 +1291 @@
 
   {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
+  {"extrapolate_second_order_sw", extrapolate_second_order_sw, METH_VARARGS, "Print out"},
   {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
   {"gravity", gravity, METH_VARARGS, "Print out"},
@@ -870 +1298 @@
    METH_VARARGS, "Print out"},
   {"h_limiter", h_limiter,
+   METH_VARARGS, "Print out"},
+  {"h_limiter_sw", h_limiter_sw,
    METH_VARARGS, "Print out"},
   {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},