Changeset 8446

Timestamp:

Jun 14, 2012, 12:31:32 AM (11 years ago)

Author:

davies

Message:

bal_dev: Updates

Location:

trunk/anuga_work/development/gareth

Files:

• experimental/balanced_dev/__init__.py (modified) (1 diff)
• experimental/balanced_dev/swb2_boundary_conditions.py (modified) (1 diff)
• experimental/balanced_dev/swb2_domain.py (modified) (1 diff)
• experimental/balanced_dev/swb2_domain_ext.c (modified) (36 diffs)
• tests/channel_floodplain/channel_floodplain1.py (modified) (4 diffs)
• tests/channel_floodplain/plotme.py (modified) (2 diffs)
• tests/dam_break/plotme.py (modified) (1 diff)
• tests/runup/runup.py (modified) (1 diff)
• tests/runup_sinusoid/runup_sinusoid.py (modified) (2 diffs)
• tests/wave/run_wave.py (modified) (2 diffs)

trunk/anuga_work/development/gareth/experimental/balanced_dev/__init__.py

@@ -9 +9 @@
 from balanced_dev.swb2_domain import Domain as Domain
 from balanced_dev.swb2_domain import Domain
-from balanced_dev.swb2_boundary_conditions import zero_mass_flux_transmissive_momentum_flux_boundary as zero_mass_flux_transmissive_momentum_flux_boundary
+#from balanced_dev.swb2_boundary_conditions import zero_mass_flux_transmissive_momentum_flux_boundary as zero_mass_flux_transmissive_momentum_flux_boundary
 #, \
 #     Transmissive_boundary,\

trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_boundary_conditions.py

@@ -6 +6 @@
 #####
 
-class  zero_mass_flux_transmissive_momentum_flux_boundary(Boundary):
-    """ Boundary which operates directly on the fluxes
-        CAN BE UNSTABLE -- CAREFUL
-        Sets boundary values of h, uh, vh as in transmissive boundary, but 
-        also ensures:
-        flux[0] = 0.0,
+#class  zero_mass_flux_transmissive_momentum_flux_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        CAN BE UNSTABLE -- CAREFUL
+#        Sets boundary values of h, uh, vh as in transmissive boundary, but 
+#        also ensures:
+#        flux[0] = 0.0,
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_transmissive_momentum_flux boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_transmissive_momentum_flux_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """Transmissive boundaries return the edge values
+#        of the volume they serve.
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#        return q
+#
+######
+#
+#class  zero_mass_flux_zero_momentum_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        Sets boundary values of h=neighbour_h, uh=0, vh=0. 
+#        and ensure that:
+#        flux[0] = 0.0
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_zero_momentum_boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_zero_momentum_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """ Apply chosen values to q[1], q[2]
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
+#
+#        r = rotate(q, normal, direction = 1)
+#        #r[1] = -r[1]
+#        tdamp=0.0
+#        ndamp=0.0
+#        r[1]=r[1]*ndamp
+#        r[2]=r[2]*tdamp
+#        q = rotate(r, normal, direction = -1)
+#
+#        # q[1]=0.
+#        # q[2]=0.
+#
+#        return q
+#
+#
+#class  zero_mass_flux_zero_n_transmissive_t_momentum_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        Sets boundary values of h=neighbour_h, uh=0, vh=vh 
+#        and ensure that:
+#        flux[0] = 0.0
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_zero_n_transmissive_t_momentum_boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_zero_n_transmissive_t_momentum_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """ Apply chosen values to q[1], q[2]
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
+#
+#        r = rotate(q, normal, direction = 1)
+#        #r[1] = -r[1]
+#        tdamp=1.0
+#        ndamp=0.0
+#        r[1]=r[1]*ndamp
+#        r[2]=r[2]*tdamp
+#        q = rotate(r, normal, direction = -1)
+#
+#        # q[1]=0.
+#        # q[2]=0.
+#
+#        return q
+#
+#
+#class  zero_mass_flux_zero_t_transmissive_n_momentum_boundary(Boundary):
+#    """ Boundary which operates directly on the fluxes
+#        Sets boundary values of h=neighbour_h, uh=neighbour_uh, vh=0.
+#        and ensure that:
+#        flux[0] = 0.0
+#    """
+#    def __init__(self, domain = None):
+#        Boundary.__init__(self)
+#
+#        if domain is None:
+#            msg = 'Domain must be specified for zero_mass_flux_zero_t_transmissive_n_momentum_boundary'
+#            raise Exception, msg
+#
+#        self.domain = domain
+#
+#    def __repr__(self):
+#        return 'zero_mass_flux_zero_t_transmissive_n_momentum_boundary(%s)' %self.domain
+#
+#    def evaluate(self, vol_id, edge_id):
+#        """ Apply chosen values to q[1], q[2]
+#        """
+#
+#        if self.domain.get_centroid_transmissive_bc() :
+#            q = self.domain.get_evolved_quantities(vol_id)
+#        else:
+#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+#
+#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
+#
+#        r = rotate(q, normal, direction = 1)
+#        #r[1] = -r[1]
+#        tdamp=0.0
+#        ndamp=1.0
+#        r[1]=r[1]*ndamp
+#        r[2]=r[2]*tdamp
+#        q = rotate(r, normal, direction = -1)
+#
+#        # q[1]=0.
+#        # q[2]=0.
+#
+#        return q
+#
+#
+class Characteristic_boundary(Boundary):
     """
-    def __init__(self, domain = None):
+
+    Example:
+
+    def waveform(t):
+        return sea_level + normalized_amplitude/cosh(t-25)**2
+
+    Bf = Nudge_boundary(domain, waveform)
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+
+    def __init__(self, domain=None, function=None):
+        """ Instantiate a
+            Nudge_boundary.
+            domain is the domain containing the boundary
+            function is the function to apply
+        """
+
         Boundary.__init__(self)
 
         if domain is None:
-            msg = 'Domain must be specified for zero_mass_flux_transmissive_momentum_flux boundary'
+            msg = 'Domain must be specified for this type boundary'
             raise Exception, msg
 
+        if function is None:
+            msg = 'Function must be specified for this type boundary'
+            raise Exception, msg
+
         self.domain = domain
+        self.function = function
+
 
     def __repr__(self):
-        return 'zero_mass_flux_transmissive_momentum_flux_boundary(%s)' %self.domain
+        """ Return a representation of this instance. """
+        msg = 'Nudge_boundary'
+        msg += '(%s)' % self.domain
+        return msg
+
 
     def evaluate(self, vol_id, edge_id):
-        """Transmissive boundaries return the edge values
-        of the volume they serve.
         """
-
-        if self.domain.get_centroid_transmissive_bc() :
-            q = self.domain.get_evolved_quantities(vol_id)
+        """
+
+        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
+        bed = self.domain.quantities['elevation'].centroid_values[vol_id]
+        depth=max(q[0]-bed,0.0)
+        dt=self.domain.timestep
+
+        normal = self.domain.get_normal(vol_id, edge_id)
+
+
+        t = self.domain.get_time()
+
+        if hasattr(self.function, 'time'):
+            # Roll boundary over if time exceeds
+            while t > self.function.time[-1]:
+                msg = 'WARNING: domain time %.2f has exceeded' % t
+                msg += 'time provided in '
+                msg += 'Flather_boundary object.\n'
+                msg += 'I will continue, reusing the object from t==0'
+                log.critical(msg)
+                t -= self.function.time[-1]
+
+        value = self.function(t)
+        try:
+            x = float(value)
+        except:
+            x = float(value[0])
+
+        ##
+        ndotq = (normal[0]*q[1] + normal[1]*q[2]) # momentum perpendicular to the boundary
+
+        if(depth==0.):
+            q[0] = x
+            q[1] = 0. 
+            q[2] = 0. 
+
         else:
-            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
+
+            # Asssume sub-critical flow. Set the values of the characteristics as
+            # appropriate, depending on whether we have inflow or outflow
+            tmp = (9.8/depth)**0.5
+            if(ndotq>0.):
+                # Only need to impose one characteristic
+                # Assume zero external velocity (dangerous?)
+                # FIXME: Allow external velocities to be specified in the function
+                # This would be needed to use the boundary for real inflows (e.g. run_wave.py)
+                # At the moment, it seems to work okay as a 'passive' outflow boundary (e.g. run_wave.py)
+                # Here it attenuates the wave at the outflow a bit, but doesn't do too much damage compared
+                # with the other boundary conditions
+                w1 = 0. - tmp*x
+                w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
+                w3 = ndotq/depth + tmp*q[0]
+                
+            else:
+                # Need to set 2 characteristics
+                # Assume zero external velocity (dangerous?)
+                # FIXME: Allow this to be specified with a function
+                w1 = 0. - tmp*x
+                w2 = 0.
+                #w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
+                w3 = ndotq/depth + tmp*q[0]
+
+
+            q[0] = (w3-w1)/(2*tmp)
+            qperp= (w3+w1)/2. *depth
+            qpar=  w2*depth
+
+            # So q[1], q[2] = qperp*(normal[0], normal[1]) + qpar*(-normal[1], normal[0])
+
+            q[1] = qperp*normal[0] + qpar*normal[1]
+            q[2] = qperp*normal[1] -qpar*normal[0]
+
         return q
-
-#####
-
-class  zero_mass_flux_zero_momentum_boundary(Boundary):
-    """ Boundary which operates directly on the fluxes
-        Sets boundary values of h=neighbour_h, uh=0, vh=0. 
-        and ensure that:
-        flux[0] = 0.0
-    """
-    def __init__(self, domain = None):
-        Boundary.__init__(self)
-
-        if domain is None:
-            msg = 'Domain must be specified for zero_mass_flux_zero_momentum_boundary'
-            raise Exception, msg
-
-        self.domain = domain
-
-    def __repr__(self):
-        return 'zero_mass_flux_zero_momentum_boundary(%s)' %self.domain
-
-    def evaluate(self, vol_id, edge_id):
-        """ Apply chosen values to q[1], q[2]
-        """
-
-        if self.domain.get_centroid_transmissive_bc() :
-            q = self.domain.get_evolved_quantities(vol_id)
-        else:
-            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
-
-        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
-
-        r = rotate(q, normal, direction = 1)
-        #r[1] = -r[1]
-        tdamp=0.0
-        ndamp=0.0
-        r[1]=r[1]*ndamp
-        r[2]=r[2]*tdamp
-        q = rotate(r, normal, direction = -1)
-
-        # q[1]=0.
-        # q[2]=0.
-
-        return q
-
-
-class  zero_mass_flux_zero_n_transmissive_t_momentum_boundary(Boundary):
-    """ Boundary which operates directly on the fluxes
-        Sets boundary values of h=neighbour_h, uh=0, vh=vh 
-        and ensure that:
-        flux[0] = 0.0
-    """
-    def __init__(self, domain = None):
-        Boundary.__init__(self)
-
-        if domain is None:
-            msg = 'Domain must be specified for zero_mass_flux_zero_n_transmissive_t_momentum_boundary'
-            raise Exception, msg
-
-        self.domain = domain
-
-    def __repr__(self):
-        return 'zero_mass_flux_zero_n_transmissive_t_momentum_boundary(%s)' %self.domain
-
-    def evaluate(self, vol_id, edge_id):
-        """ Apply chosen values to q[1], q[2]
-        """
-
-        if self.domain.get_centroid_transmissive_bc() :
-            q = self.domain.get_evolved_quantities(vol_id)
-        else:
-            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
-
-        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
-
-        r = rotate(q, normal, direction = 1)
-        #r[1] = -r[1]
-        tdamp=1.0
-        ndamp=0.0
-        r[1]=r[1]*ndamp
-        r[2]=r[2]*tdamp
-        q = rotate(r, normal, direction = -1)
-
-        # q[1]=0.
-        # q[2]=0.
-
-        return q
-
-
-class  zero_mass_flux_zero_t_transmissive_n_momentum_boundary(Boundary):
-    """ Boundary which operates directly on the fluxes
-        Sets boundary values of h=neighbour_h, uh=neighbour_uh, vh=0.
-        and ensure that:
-        flux[0] = 0.0
-    """
-    def __init__(self, domain = None):
-        Boundary.__init__(self)
-
-        if domain is None:
-            msg = 'Domain must be specified for zero_mass_flux_zero_t_transmissive_n_momentum_boundary'
-            raise Exception, msg
-
-        self.domain = domain
-
-    def __repr__(self):
-        return 'zero_mass_flux_zero_t_transmissive_n_momentum_boundary(%s)' %self.domain
-
-    def evaluate(self, vol_id, edge_id):
-        """ Apply chosen values to q[1], q[2]
-        """
-
-        if self.domain.get_centroid_transmissive_bc() :
-            q = self.domain.get_evolved_quantities(vol_id)
-        else:
-            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
-
-        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
-
-        r = rotate(q, normal, direction = 1)
-        #r[1] = -r[1]
-        tdamp=0.0
-        ndamp=1.0
-        r[1]=r[1]*ndamp
-        r[2]=r[2]*tdamp
-        q = rotate(r, normal, direction = -1)
-
-        # q[1]=0.
-        # q[2]=0.
-
-        return q
-
-
 ###
 def rotate(q,normal,direction=1):

trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py

@@ -75 +75 @@
         # Note that the extrapolation method used in quantity_ext.c (e.g.
         # extrapolate_second_order_and_limit_by_edge) uses a constant value for
-        # all the betas == beta_rk2 in our case.  
-        self.beta_w=2.0
-        self.beta_w_dry=2.0
-        self.beta_uh=2.0
-        self.beta_uh_dry=2.0
-        self.beta_vh=2.0
-        self.beta_vh_dry=2.0
+        # all the betas.  
+        self.beta_w=1.0
+        self.beta_w_dry=1.0
+        self.beta_uh=1.0
+        self.beta_uh_dry=1.0
+        self.beta_vh=1.0
+        self.beta_vh_dry=1.0
 
         #self.optimise_dry_cells=True

trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

@@ -227 +227 @@
       edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
       //if(i<2) {
-        edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
+      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
       //}else{
       //  edgeflux[i] += 0.5*s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
@@ -307 +307 @@
     int useint;
     double stage_edge_lim, scale_factor_shallow, bedtop, bedbot, pressure_flux, hc, hc_n;
-    double max_bed_edgevalue[number_of_elements];
-    //double depthtemp, max_bed_edge_value; 
+    //double min_bed_edgevalue[number_of_elements]; // FIXME: Set memory explicitly
     static long call = 1; // Static local variable flagging already computed flux
 
+    double *min_bed_edgevalue;
+
+    min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
     // Start computation
     call++; // Flag 'id' of flux calculation for this timestep
@@ -321 +323 @@
 
 
-    // Compute maximum bed edge value on each triangle
-    //for (k = 0; k < number_of_elements; k++){
-    //    max_bed_edgevalue[k] = max(bed_edge_values[3*k], 
-    //                               max(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
-    //
-    //}
+    // Compute minimum bed edge value on each triangle
+    for (k = 0; k < number_of_elements; k++){
+        min_bed_edgevalue[k] = min(bed_edge_values[3*k], 
+                                   min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
+        //min_bed_edgevalue[k] = bed_centroid_values[k];
+    
+    }
 
 
@@ -402 +405 @@
             // NOTE: This also means no bedslope term -- which is arguably
             // appropriate, since the depth is zero at the cell centroid
-            if(( hc == 0.0)&((hc_n == 0.0)|((n<0)&(qr[0]<zl)))){
+            if(( hc == 0.0)&(hc_n == 0.0)){
+                already_computed_flux[ki] = call; // #k Done
+                already_computed_flux[nm] = call; // #n Done
+                max_speed = 0.0;
                 continue;
             }
@@ -410 +416 @@
             if(n>=0){
                 if(hc==0.0){
-                    ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
-                    //ql[0]=qr[0]; // max(min(qr[0],stage_centroid_values[k]),zl);
+                    //ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
+                    //ql[0]=max(min(qr[0],0.5*(stage_centroid_values[k]+stage_centroid_values[n])),zl);
+                    ql[0]=qr[0]; // max(min(qr[0],stage_centroid_values[k]),zl);
                 }
                 if(hc_n==0.0){
-                    qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
-                    //qr[0]=ql[0]; 
+                    //qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
+                    //qr[0]=max(min(ql[0],0.5*(stage_centroid_values[n]+stage_centroid_values[k])),zr);
+                    qr[0]=ql[0]; 
                 }
             }else{
@@ -421 +429 @@
                 if((hc==0.0)){
                     ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
+                    //ql[0]=max(min(qr[0],ql[0]),zl);
                 }
             }
@@ -432 +441 @@
                     epsilon, h0, limiting_threshold, g,
                     edgeflux, &max_speed, &pressure_flux);
-            //_flux_function_fraccaro(ql, qr, zl, zr,
-            //        normals[ki2], normals[ki2 + 1],
-            //        epsilon, h0, limiting_threshold, g,
-            //        edgeflux, &max_speed);
-   
+
+            // Prevent outflow from 'seriously' dry cells
+            if(stage_centroid_values[k]<min_bed_edgevalue[k]){
+                if(edgeflux[0]>0.0){
+                    edgeflux[0]=0.0;
+                    edgeflux[1]=0.0;
+                    edgeflux[2]=0.0;
+                }
+            }
+            if(n>=0){
+                if(stage_centroid_values[n]<min_bed_edgevalue[n]){
+                    if(edgeflux[0]<0.0){
+                        edgeflux[0]=0.0;
+                        edgeflux[1]=0.0;
+                        edgeflux[2]=0.0;
+                    }
+                }
+            } 
+
             // Multiply edgeflux by edgelength
             length = edgelengths[ki];
@@ -443 +466 @@
             edgeflux[2] *= length;
 
-            if(n<0){
-                if(boundary_flux_type[m]>0){
-                    // IMPOSE BOUNDARY CONDITIONS DIRECTLY ON THE FLUX
-                    if(boundary_flux_type[m]==1){
-                        // Zero_mass_flux_transmissive_momentum_flux boundary, or
-                        // zero_mass_flux_zero_momentum_boundary
-                        // Just need to set the mass flux to zero, as the
-                        // transmissive momentum is ensured by the values of
-                        // qr[0], qr[1], qr[2]
-                        edgeflux[0] = 0.0;
-                    }
-                }
-            }
+            // GET RID OF THIS: EXPERIMENTAL
+            //if(n<0){
+            //    if(boundary_flux_type[m]>0){
+            //        // IMPOSE BOUNDARY CONDITIONS DIRECTLY ON THE FLUX
+            //        if(boundary_flux_type[m]==1){
+            //            // Zero_mass_flux_transmissive_momentum_flux boundary, or
+            //            // zero_mass_flux_zero_momentum_boundary
+            //            // Just need to set the mass flux to zero, as the
+            //            // transmissive momentum is ensured by the values of
+            //            // qr[0], qr[1], qr[2]
+            //            printf("Warning: WEIRD EDGEFLUX THING IS ACTING");
+            //            edgeflux[0] = 0.0;
+            //        }
+            //    }
+            //}
 
             // Update triangle k with flux from edge i
@@ -462 +487 @@
             ymom_explicit_update[k] -= edgeflux[2];
 
-            // Calculate the bed slope term, using the 'divergence form' from 
-            // Alessandro Valiani and Lorenzo Begnudelli, Divergence Form for Bed Slope Source Term in Shallow Water Equations
-            // Journal of Hydraulic Engineering, 2006, 132:652-665
-            if(hc>0.0e-03){
+            // Compute bed slope term
+            if(hc>-9999.0){
                 //Bedslope approx 1:
                 //bedslope_work = g*hc*(ql[0])*length-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length;
@@ -486 +509 @@
             }else{
                 // Treat nearly dry cells 
-                //bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
+                bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
                 //
-                bedslope_work = -pressure_flux*length;
+                //bedslope_work = -pressure_flux*length;
                 xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
                 ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
@@ -506 +529 @@
                 ymom_explicit_update[n] += edgeflux[2];
                 //Add bed slope term here
-                if(hc_n>0.0e-03){
+                if(hc_n>-9999.0){
                 //if(stage_centroid_values[n] > max_bed_edgevalue[n]){
                     // Bedslope approx 1:
@@ -528 +551 @@
                 }else{
                     // Treat nearly dry cells
-                    //bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
-                    bedslope_work = -pressure_flux*length;
+                    bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
+                    //bedslope_work = -pressure_flux*length;
                     xmom_explicit_update[n] += normals[ki2]*bedslope_work;
                     ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
@@ -578 +601 @@
     } // End triangle k
 
+    free(min_bed_edgevalue);
 
     return timestep;
@@ -599 +623 @@
 
   // This acts like minimum_allowed height, but scales with the vertical
-  // distance between the bed_centroid_value and the max bed_edge_value of every triangle.
-  double minimum_relative_height=0.05; 
+  // distance between the bed_centroid_value and the max bed_edge_value of
+  // every triangle.
+  double minimum_relative_height=0.1; 
 
 
@@ -619 +644 @@
         if (hc <= 0.0){
              // Definine the minimum bed edge value on triangle k.
-             bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
+             // Setting = minimum edge value can lead to mass conservation problems
+             //bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
+             //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
+             // Setting = minimum vertex value seems better, but might tend to be less smooth 
+             bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
              // Minimum allowed stage = bmin
+             // WARNING: ADDING MASS if wc[k]<bmin
+             if(wc[k]<bmin){
+                 printf("Adding mass to dry cell\n");
+             }
              wc[k] = max(wc[k], bmin) ; 
+
+             // Set vertex values as well. Seems that this shouldn't be
+             // needed. However from memory this is important at the first
+             // time step, for 'dry' areas where the designated stage is
+             // less than the bed centroid value
              wv[3*k] = max(wv[3*k], bmin);
              wv[3*k+1] = max(wv[3*k+1], bmin);
@@ -683 +721 @@
   // dq2=q(vertex2)-q(vertex0)
 
-  // This is a simple implementation -- the one below is the
-  // original version -- I wonder if it is any faster?
+  // This is a simple implementation 
   *qmax = max(max(dq0, max(dq0+dq1, dq0+dq2)), 0.0) ;
   *qmin = min(min(dq0, min(dq0+dq1, dq0+dq2)), 0.0) ;
  
-  //if (dq0>=0.0){
-  //  if (dq1>=dq2){
-  //    if (dq1>=0.0)
-  //  *qmax=dq0+dq1;
-  //    else
-  //  *qmax=dq0;
-  //    
-  //    *qmin=dq0+dq2;
-  //    if (*qmin>=0.0) *qmin = 0.0;
-  //  }
-  //  else{// dq1<dq2
-  //    if (dq2>0)
-  //  *qmax=dq0+dq2;
-  //    else
-  //  *qmax=dq0;
-  //    
-  //    *qmin=dq0+dq1;    
-  //    if (*qmin>=0.0) *qmin=0.0;
-  //  }
-  //}
-  //else{//dq0<0
-  //  if (dq1<=dq2){
-  //    if (dq1<0.0)
-  //  *qmin=dq0+dq1;
-  //    else
-  //  *qmin=dq0;
-  //    
-  //    *qmax=dq0+dq2;    
-  //    if (*qmax<=0.0) *qmax=0.0;
-  //  }
-  //  else{// dq1>dq2
-  //    if (dq2<0.0)
-  //  *qmin=dq0+dq2;
-  //    else
-  //  *qmin=dq0;
-  //    
-  //    *qmax=dq0+dq1;
-  //    if (*qmax<=0.0) *qmax=0.0;
-  //  }
-  //}
   return 0;
 }
@@ -759 +756 @@
   phi=min(r*beta_w,1.0);
   //phi=1.;
-  //for (i=0;i<3;i++)
   dqv[0]=dqv[0]*phi;
   dqv[1]=dqv[1]*phi;
@@ -766 +762 @@
   return 0;
 }
+
+//int limit_gradient2(double *dqv, double qmin, double qmax, double beta_w, double r0scale){
+//  // EXPERIMENTAL LIMITER, DOESN'T WORK 
+//  // 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.
+//  // FIXME: Perhaps use the epsilon used elsewhere.
+//  
+//  for (i=0;i<3;i++){
+//    if (dqv[i]<-TINY)
+//      r0=qmin/dqv[i]*r0scale*2.0;
+//      
+//    if (dqv[i]>TINY)
+//      r0=qmax/dqv[i]*r0scale*2.0;
+//      
+//    r=min(r0,r);
+//  }
+//  
+//  phi=min(r*beta_w,1.0);
+//  //phi=1.;
+//  //for (i=0;i<3;i++)
+//  dqv[0]=dqv[0]*phi;
+//  dqv[1]=dqv[1]*phi;
+//  dqv[2]=dqv[2]*phi;
+//    
+//  return 0;
+//}
 
 // Computational routine
@@ -789 +819 @@
                                  double* ymom_edge_values,
                                  double* elevation_edge_values,
+                                 double* stage_vertex_values,
+                                 double* xmom_vertex_values,
+                                 double* ymom_vertex_values,
+                                 double* elevation_vertex_values,
                                  int optimise_dry_cells, 
                                  int extrapolate_velocity_second_order) {
                   
-  // Note that although this routine refers to EDGE values, it can equally well
-  // be applied to vertex values -- it is a modified version of such a routine.
-  // Also note that momentum can be extrapolated using velocity, via
-  // modifications at the start and end of the routine.                 
-
   // Local variables
   double a, b; // Gradient vector used to calculate edge values from centroids
@@ -804 +833 @@
   double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
   double hc, h0, h1, h2, beta_tmp, hfactor;
-  double dk, dv0, dv1, dv2, de[3];
-  //double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements]
-  //double stage_centroid_store[number_of_elements];
+  double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale;
   
   double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store;
@@ -828 +855 @@
           ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
 
-          // Experimental -- replace stage with a weighting of 'stage' and
-          // 'depth', where the weights depend on the froude number **2
-          // stage_centroid_store[k]=stage_centroid_values[k];
-          // hfactor =(xmom_centroid_values[k]*xmom_centroid_values[k] + 
-          //          ymom_centroid_values[k]*ymom_centroid_values[k])/20.0;
-          //          //(9.8*max(stage_centroid_values[k] - elevation_centroid_values[k],minimum_allowed_height));
-          //stage_centroid_values[k] -= min(1.0, hfactor)*elevation_centroid_values[k];
       }
   }
@@ -887 +907 @@
       dyv1 = yv1 - y; 
       dyv2 = yv2 - y;
+      // Compute the minimum distance from the centroid to an edge
+      //demin=min(dxv0*dxv0 +dyv0*dyv0, min(dxv1*dxv1+dyv1*dyv1, dxv2*dxv2+dyv2*dyv2));
+      //demin=sqrt(demin);
     }
 
@@ -911 +934 @@
       // triangle area will be zero, and a first order extrapolation will be
       // used
-      //if(stage_centroid_values[k2]<=elevation_centroid_values[k2]){
-      //    k2 = k ;
-      //}
-      //if(stage_centroid_values[k0]<=elevation_centroid_values[k0]){
-      //    k0 = k ;
-      //}
-      //if(stage_centroid_values[k1]<=elevation_centroid_values[k1]){
-      //    k1 = k ;
-      //}
+      if(stage_centroid_values[k2]<=elevation_centroid_values[k2]){
+          k2 = k ;
+      }
+      if(stage_centroid_values[k0]<=elevation_centroid_values[k0]){
+          k0 = k ;
+      }
+      if(stage_centroid_values[k1]<=elevation_centroid_values[k1]){
+          k1 = k ;
+      }
       // Alternative approach (BRUTAL) -- if the max neighbour bed elevation is greater
       // than the min neighbour stage, then we use first order extrapolation
-      bedmax = max(elevation_centroid_values[k], 
-                   max(elevation_centroid_values[k0],
-                       max(elevation_centroid_values[k1], elevation_centroid_values[k2])));
-      stagemin =min(stage_centroid_values[k], 
-                   min(stage_centroid_values[k0],
-                       min(stage_centroid_values[k1], stage_centroid_values[k2])));
-  
-      if(stagemin < bedmax){
-         // This will cause first order extrapolation
-         k2 = k;
-         k0 = k;
-         k1 = k;
-      } 
+      //bedmax = max(elevation_centroid_values[k], 
+      //             max(elevation_centroid_values[k0],
+      //                 max(elevation_centroid_values[k1], elevation_centroid_values[k2])));
+      //stagemin = min(stage_centroid_values[k], 
+      //             min(stage_centroid_values[k0],
+      //                 min(stage_centroid_values[k1], stage_centroid_values[k2])));
+      // 
+      //if(stagemin < bedmax){
+      //   // This will cause first order extrapolation
+      //   k2 = k;
+      //   k0 = k;
+      //   k1 = k;
+      //} 
      
       // Get the auxiliary triangle's vertex coordinates 
@@ -949 +972 @@
       x2 = centroid_coordinates[coord_index];
       y2 = centroid_coordinates[coord_index+1];
-      
+
+      // compute the maximum distance from the centroid to a neighbouring
+      // centroid
+      //dcmax=max( (x0-x)*(x0-x) + (y0-y)*(y0-y),     
+      //           max((x1-x)*(x1-x) + (y1-y)*(y1-y), 
+      //               (x2-x)*(x2-x) + (y2-y)*(y2-y)));
+      //dcmax=sqrt(dcmax);
+      //// Ratio of centroid to edge distance -- useful in attempting to adapt limiter
+      //if(dcmax>0.){
+      //    r0scale=demin/dcmax; 
+      //    //printf("%f \n", r0scale);
+      //}else{
+      //    r0scale=0.5;
+      //}
+
       // Store x- and y- differentials for the vertices 
       // of the auxiliary triangle
@@ -970 +1007 @@
           //return -1;
 
-          //stage_edge_values[k3]   = 0.5*(stage_centroid_values[k]+stage_centroid_values[k0]);
-          //stage_edge_values[k3+1] = 0.5*(stage_centroid_values[k]+stage_centroid_values[k1]);
-          //stage_edge_values[k3+2] = 0.5*(stage_centroid_values[k]+stage_centroid_values[k2]);
           stage_edge_values[k3]   = stage_centroid_values[k];
           stage_edge_values[k3+1] = stage_centroid_values[k];
@@ -1002 +1036 @@
       {
         hfactor = 1.0 ;//hmin/(hmin + 0.004);
+        //hfactor=hmin/(hmin + 0.004);
       }
       
@@ -1048 +1083 @@
       find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
       
-      // Playing with dry wet interface
-      //hmin = qmin;
-      //beta_tmp = beta_w_dry;
-      //if (hmin>minimum_allowed_height)
       beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
     
-      //printf("min_alled_height = %f\n",minimum_allowed_height);
-      //printf("hmin = %f\n",hmin);
-      //printf("beta_w = %f\n",beta_w);
-      //printf("beta_tmp = %f\n",beta_tmp);
       // Limit the gradient
       limit_gradient(dqv, qmin, qmax, beta_tmp); 
+      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale); 
       
       //for (i=0;i<3;i++)
@@ -1066 +1094 @@
       stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
       
-      //if(k==300){
-      //  printf("Stage centroid values: %f, v1, %f, v2, %f, v3, %f \n", stage_centroid_values[k], 
-      //          stage_edge_values[k3+0], stage_edge_values[k3+1], stage_edge_values[k3+2]);
-      //}     
- 
       //-----------------------------------
       // xmomentum
@@ -1102 +1125 @@
       find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
 
-      //beta_tmp = beta_uh;
-      //if (hmin<minimum_allowed_height)
-      //beta_tmp = beta_uh_dry;
       beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
-
-      //if(k==300){
-      //  printf("dqv0: %f, dqv1, %f, dqv2, %f \n", dqv[0],dqv[1], dqv[2]) ;
-
-      //}
 
       // Limit the gradient
       limit_gradient(dqv, qmin, qmax, beta_tmp);
-      
-      //if(k==300){
-      //  printf("dqv0: %f, dqv1, %f, dqv2, %f \n", dqv[0],dqv[1], dqv[2]) ;
-
-      //}
+      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
+      
 
       for (i=0; i < 3; i++)
@@ -1124 +1136 @@
         xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
       }
-      
-      //if(k==300){
-      //  printf("xmom centroid values: %f, v1, %f, v2, %f, v3, %f \n", xmom_centroid_values[k], 
-      //          xmom_edge_values[k3+0], xmom_edge_values[k3+1], xmom_edge_values[k3+2]);
-      //}     
       
       //-----------------------------------
@@ -1161 +1168 @@
       find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
       
-      //beta_tmp = beta_vh;
-      //
-      //if (hmin<minimum_allowed_height)
-      //beta_tmp = beta_vh_dry;
       beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
 
       // Limit the gradient
       limit_gradient(dqv, qmin, qmax, beta_tmp);
+      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
       
       for (i=0;i<3;i++)
@@ -1175 +1179 @@
       }
       
-      //if(k==300){
-      //  printf("ymom centroid values: %f, v1, %f, v2, %f, v3, %f \n", ymom_centroid_values[k], 
-      //          ymom_edge_values[k3+0], ymom_edge_values[k3+1], ymom_edge_values[k3+2]);
-      //}     
     } // End number_of_boundaries <=1 
     else
@@ -1346 +1346 @@
       limit_gradient(dqv, qmin, qmax, beta_w);
       
-      //for (i=0;i<3;i++)
-      //ymom_edge_values[k3] = ymom_centroid_values[k] + dqv[0];
-      //ymom_edge_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
-      //ymom_edge_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
-
-for (i=0;i<3;i++)
-        {
-        ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
-        }
-//ymom_edge_values[k3] = ymom_centroid_values[k] + dqv[0];
-//ymom_edge_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
-//ymom_edge_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
+      for (i=0;i<3;i++)
+              {
+              ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
+              }
     } // else [number_of_boundaries==2]
   } // for k=0 to number_of_elements-1
   
-  if(extrapolate_velocity_second_order==1){
-      // Convert back from velocity to momentum
-      
-      for (k=0; k<number_of_elements; k++){
-          k3=3*k;
-
-          // Convert energy back to stage
-          //stage_centroid_values[k] = stage_centroid_store[k];
-
-          //Convert velocity back to momenta
+  // Compute vertex values of quantities
+  for (k=0; k<number_of_elements; k++){
+      k3=3*k;
+
+      // Compute stage vertex values 
+      stage_vertex_values[k3] = stage_edge_values[k3+1] + stage_edge_values[k3+2] -stage_edge_values[k3] ; 
+      stage_vertex_values[k3+1] =  stage_edge_values[k3] + stage_edge_values[k3+2]-stage_edge_values[k3+1]; 
+      stage_vertex_values[k3+2] =  stage_edge_values[k3] + stage_edge_values[k3+1]-stage_edge_values[k3+2]; 
+      
+     // Compute xmom vertex values 
+      xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ; 
+      xmom_vertex_values[k3+1] =  xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1]; 
+      xmom_vertex_values[k3+2] =  xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2]; 
+
+      // Compute ymom vertex values 
+      ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ; 
+      ymom_vertex_values[k3+1] =  ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1]; 
+      ymom_vertex_values[k3+2] =  ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2]; 
+
+      // If needed, convert from velocity to momenta
+      if(extrapolate_velocity_second_order==1){
+          //Convert velocity back to momenta at centroids
           xmom_centroid_values[k] = xmom_centroid_store[k];
           ymom_centroid_values[k] = ymom_centroid_store[k];
 
+          // Re-compute momenta at edges
           for (i=0; i<3; i++){
-              //stage_edge_values[k3+i] -=(xmom_edge_values[k3+i]*xmom_edge_values[k3+i] 
-              //                           + ymom_edge_values[k3+i]*ymom_edge_values[k3+i])*0.0;
-              //stage_edge_values[k3+i] += elevation_edge_values[k3+i];
-              //de[i] = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);
-              //for(ii=0; ii<1; ii++){
-              //hfactor =(xmom_edge_values[k3+i]*xmom_edge_values[k3+i] + 
-              //          ymom_edge_values[k3+i]*ymom_edge_values[k3+i])/20.0; // 20 ~= 2*g
-                        //(9.8*de[i]);
-              //hfactor = 0.0; //max(hfactor/20.0 - 0.2, 0.0);
-              //stage_edge_values[k3+i] = stage_edge_values[k3+i];//+ min(hfactor, 1.0)*elevation_edge_values[k3+i];
-              //de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],1.0e-03);
-              //}
-
               de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0);
               xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i];
@@ -1393 +1385 @@
           }
 
-          
-      } 
-  }
+          // Re-compute momenta at vertices
+          for (i=0; i<3; i++){
+              de[i] = max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0);
+              xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i];
+              ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i];
+          }
+
+      }
+
+   
+  } 
 
   free(xmom_centroid_store);
@@ -1404 +1404 @@
 }           
 
-int extrapolate_from_edges_to_vertices(int number_of_elements,
-                                        double* stage_edge_values,
-                                        double* xmom_edge_values,
-                                        double* ymom_edge_values,
-                                        double* elevation_edge_values,
-                                        double* stage_vertex_values,
-                                        double* xmom_vertex_values,
-                                        double* ymom_vertex_values,
-                                        double* elevation_vertex_values,
-                                        int extrapolate_velocity_second_order) {
-
-  int i, i3;
-  double v0, v1, v2;
-
-  for(i=0; i<number_of_elements; i++){
-        i3 = 3*i;
-        stage_vertex_values[i3] = stage_edge_values[i3+1] + stage_edge_values[i3+2] -stage_edge_values[i3] ; 
-        stage_vertex_values[i3+1] =  stage_edge_values[i3] + stage_edge_values[i3+2]-stage_edge_values[i3+1]; 
-        stage_vertex_values[i3+2] =  stage_edge_values[i3] + stage_edge_values[i3+1]-stage_edge_values[i3+2]; 
-
-        if(extrapolate_velocity_second_order==1){
-            // Compute x-velocity, carefully to avoid division by zero
-            //v0 = xmom_edge_values[i3+0]/(stage_edge_values[i3+0]-elevation_edge_values[i3+0]); 
-            //v1 = xmom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]); 
-            //v2 = xmom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
-            if(stage_edge_values[i3]>elevation_edge_values[i3]){
-                v0 = xmom_edge_values[i3]/(stage_edge_values[i3]-elevation_edge_values[i3]);
-            }else{
-                v0 = 0.0;
-            }
-            if(stage_edge_values[i3+1]>elevation_edge_values[i3+1]){
-                v1 = xmom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);
-            }else{
-                v1 = 0.0;
-            }
-
-            if(stage_edge_values[i3+2]>elevation_edge_values[i3+2]){
-                v2 = xmom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
-            }else{
-                v2 = 0.0;
-            }
-            //v1 = xmom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]); 
-            //v2 = xmom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
-
-            xmom_vertex_values[i3] = (-v0 + v1 +v2)*(stage_vertex_values[i3]-elevation_vertex_values[i3]);
-            xmom_vertex_values[i3+1] = (-v1 + v0 +v2)*(stage_vertex_values[i3+1]-elevation_vertex_values[i3+1]);
-            xmom_vertex_values[i3+2] = (-v2 + v0 +v1)*(stage_vertex_values[i3+2]-elevation_vertex_values[i3+2]);
-            
-            // Compute y-velocity, carefully to avoid division by zero
-            //v0 = ymom_edge_values[i3+0]/(stage_edge_values[i3+0]-elevation_edge_values[i3+0]); 
-            //v1 = ymom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]); 
-            //v2 = ymom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
-            if(stage_edge_values[i3]>elevation_edge_values[i3]){
-                v0 = ymom_edge_values[i3]/(stage_edge_values[i3]-elevation_edge_values[i3]);
-            }else{
-                v0 = 0.0;
-            }
-            if(stage_edge_values[i3+1]>elevation_edge_values[i3+1]){
-                v1 = ymom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);
-            }else{
-                v1 = 0.0;
-            }
-
-            if(stage_edge_values[i3+2]>elevation_edge_values[i3+2]){
-                v2 = ymom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
-            }else{
-                v2 = 0.0;
-            }
-            //v0 = ymom_edge_values[i3]/(stage_edge_values[i3]-elevation_edge_values[i3]); 
-            //v1 = ymom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]); 
-            //v2 = ymom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
-
-            ymom_vertex_values[i3] = (-v0 + v1 +v2)*(stage_vertex_values[i3]-elevation_vertex_values[i3]);
-            ymom_vertex_values[i3+1] = (-v1 + v0 +v2)*(stage_vertex_values[i3+1]-elevation_vertex_values[i3+1]);
-            ymom_vertex_values[i3+2] = (-v2 + v0 +v1)*(stage_vertex_values[i3+2]-elevation_vertex_values[i3+2]);
-
-        }else{
-            
-            xmom_vertex_values[i3] = xmom_edge_values[i3+1] + xmom_edge_values[i3+2]-xmom_edge_values[i3] ; 
-            xmom_vertex_values[i3+1] =  xmom_edge_values[i3] + xmom_edge_values[i3+2]-xmom_edge_values[i3+1]; 
-            xmom_vertex_values[i3+2] =  xmom_edge_values[i3] + xmom_edge_values[i3+1]-xmom_edge_values[i3+2]; 
-
-            ymom_vertex_values[i3] =  ymom_edge_values[i3+1] + ymom_edge_values[i3+2]-ymom_edge_values[i3]; 
-            ymom_vertex_values[i3+1] =  ymom_edge_values[i3] + ymom_edge_values[i3+2] -ymom_edge_values[i3+1]; 
-            ymom_vertex_values[i3+2] =  ymom_edge_values[i3] + ymom_edge_values[i3+1]-ymom_edge_values[i3+2] ; 
-
-        }
-
-    }
-    return 0 ;
-}
 //=========================================================================
 // Python Glue
 //=========================================================================
 
-
-//PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
-//  //
-//  // r = rotate(q, normal, direction=1)
-//  //
-//  // Where q is assumed to be a Float numeric array of length 3 and
-//  // normal a Float numeric array of length 2.
-//
-//  // FIXME(Ole): I don't think this is used anymore
-//
-//  PyObject *Q, *Normal;
-//  PyArrayObject *q, *r, *normal;
-//
-//  static char *argnames[] = {"q", "normal", "direction", NULL};
-//  int dimensions[1], i, direction=1;
-//  double n1, n2;
-//
-//  // Convert Python arguments to C
-//  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
-//                   &Q, &Normal, &direction)) {
-//    report_python_error(AT, "could not parse input arguments");
-//    return NULL;
-//  }  
-//
-//  // Input checks (convert sequences into numeric arrays)
-//  q = (PyArrayObject *)
-//    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
-//  normal = (PyArrayObject *)
-//    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
-//
-//
-//  if (normal -> dimensions[0] != 2) {
-//    report_python_error(AT, "Normal vector must have 2 components");
-//    return NULL;
-//  }
-//
-//  // Allocate space for return vector r (don't DECREF)
-//  dimensions[0] = 3;
-//  r = (PyArrayObject *) anuga_FromDims(1, dimensions, PyArray_DOUBLE);
-//
-//  // Copy
-//  for (i=0; i<3; i++) {
-//    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
-//  }
-//
-//  // Get normal and direction
-//  n1 = ((double *) normal -> data)[0];
-//  n2 = ((double *) normal -> data)[1];
-//  if (direction == -1) n2 = -n2;
-//
-//  // Rotate
-//  _rotate((double *) r -> data, n1, n2);
-//
-//  // Release numeric arrays
-//  Py_DECREF(q);
-//  Py_DECREF(normal);
-//
-//  // Return result using PyArray to avoid memory leak
-//  return PyArray_Return(r);
-//}
 
 //========================================================================
@@ -1981 +1830 @@
                    (double*) ymom_edge_values -> data,
                    (double*) elevation_edge_values -> data,
+                   (double*) stage_vertex_values -> data,
+                   (double*) xmom_vertex_values -> data,
+                   (double*) ymom_vertex_values -> data,
+                   (double*) elevation_vertex_values -> data,
                    optimise_dry_cells, 
                    extrapolate_velocity_second_order);
-
-  //printf("In C before edges-to-vertices");
-
-  //Extrapolate from edges to vertices
-  e2 = extrapolate_from_edges_to_vertices(number_of_elements,
-                                      (double *) stage_edge_values -> data,
-                                      (double *) xmom_edge_values -> data,
-                                      (double *) ymom_edge_values -> data,
-                                      (double *) elevation_edge_values -> data,
-                                      (double *) stage_vertex_values -> data,
-                                      (double *) xmom_vertex_values -> data,
-                                      (double *) ymom_vertex_values -> data,
-                                      (double *) elevation_vertex_values -> data, 
-                                      extrapolate_velocity_second_order);
-  //printf("In C after edges-to-vertices");
 
   if (e == -1) {

trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py

@@ -10 +10 @@
 import numpy
 from anuga.structures.inlet_operator import Inlet_operator
-from anuga import *
+#from anuga import *
 #from swb_domain import domain
 #from anuga import *
@@ -82 +82 @@
 
 
-domain.set_name('channel_floodplain1_bal_dev_lowvisc') # Output name
+domain.set_name('channel_floodplain1') # Output name
 #domain.set_store_vertices_uniquely(True)
 #domain.use_edge_limiter=False
@@ -142 +142 @@
 
 # Define inlet operator 
-flow_in_yval=100.0
+flow_in_yval=10.0
 if True:
     line1 = [ [floodplain_width/2. - chan_width/2., flow_in_yval],\
@@ -245 +245 @@
     if( numpy.floor(t/100.) == t/100. ):
         print '#### COMPUTING FLOW THROUGH CROSS-SECTIONS########'
+        s0 = domain.get_flow_through_cross_section([[0., 10.0], [floodplain_width, 10.]])
         s1 = domain.get_flow_through_cross_section([[0., floodplain_length-300.0], [floodplain_width, floodplain_length-300.0]])
         s2 = domain.get_flow_through_cross_section([[0., floodplain_length-1.0], [floodplain_width, floodplain_length-1.0]])
         
-        print 'Cross sectional flows: ', s1, s2 
+        print 'Cross sectional flows: ', s0, s1, s2 
         print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
 

trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py

@@ -4 +4 @@
 
 # Time-index to plot outputs from
-index=53
+index=60
 
 #p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)
@@ -10 +10 @@
 #p=util.get_centroids(p2, velocity_extrapolation=True)
 
-p2 = util.get_output('channel_floodplain1_bal_dev_lowvisc.sww', 0.01)
+#p2 = util.get_output('channel_floodplain1_bal_dev_lowvisc.sww', 0.01)
+p2 = util.get_output('channel_floodplain1.sww', 0.01)
 p=util.get_centroids(p2, velocity_extrapolation=True)
 

trunk/anuga_work/development/gareth/tests/dam_break/plotme.py

@@ -9 +9 @@
 p2_st=util.get_centroids(p_st)
 
-
-p_dev = util.get_output('dam_break_20120404_230353/dam_break.sww', 0.001)
+p_dev = util.get_output('dam_break_20120517_003225/dam_break.sww', 0.001)
 p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
 

trunk/anuga_work/development/gareth/tests/runup/runup.py

@@ -67 +67 @@
 #------------------------------
 
-#xwrite=open("xvel.out","wb")
-#ywrite=open("yvel.out","wb")
-## Set print options to be compatible with file writing via the 'print' statement 
-#numpy.set_printoptions(threshold=numpy.nan, linewidth=numpy.nan)
-
-for t in domain.evolve(yieldstep=0.2,finaltime=30.0):
+for t in domain.evolve(yieldstep=0.2,finaltime=60.0):
     print domain.timestepping_statistics()
-
+    uh=domain.quantities['xmomentum'].centroid_values
+    vh=domain.quantities['ymomentum'].centroid_values
+    depth=domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
+    depth=depth*(depth>1.0e-06) + 1.0e-06
+    vel=((uh/depth)**2 + (vh/depth)**2)**0.5
+    print 'peak speed is', vel.max()
 
 print 'Finished'

trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

@@ -65 +65 @@
 Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
 Bd=anuga.Dirichlet_boundary([-0.1*scale_me,0.,0.])       # Constant boundary values -- not used in this example
+def waveform(t):
+    return (0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1
+Bt2=anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain,waveform)
 #Bw=anuga.Time_boundary(domain=domain,
 #       f=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,0.0,0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.
@@ -82 +85 @@
     yy = domain.quantities['ymomentum'].centroid_values
     dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values 
-    dd = (dd)*(dd>1.0e-03)+1.0e-06
+    dd = (dd)*(dd>1.0e-06)+1.0e-03
     vv = ( (xx/dd)**2 + (yy/dd)**2 )**0.5
     vv = vv*(dd>1.0e-03)

trunk/anuga_work/development/gareth/tests/wave/run_wave.py

@@ -81 +81 @@
 Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
 Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary
-Bz = swb2_boundary_conditions.zero_mass_flux_zero_momentum_boundary(domain) # Strong reflections
-Bz1 = swb2_boundary_conditions.zero_mass_flux_zero_t_transmissive_n_momentum_boundary(domain) # Strong reflections
-Bz2 = swb2_boundary_conditions.zero_mass_flux_zero_n_transmissive_t_momentum_boundary(domain) # Strong reflections
+#Bz = swb2_boundary_conditions.zero_mass_flux_zero_momentum_boundary(domain) # Strong reflections
+#Bz1 = swb2_boundary_conditions.zero_mass_flux_zero_t_transmissive_n_momentum_boundary(domain) # Strong reflections
+#Bz2 = swb2_boundary_conditions.zero_mass_flux_zero_n_transmissive_t_momentum_boundary(domain) # Strong reflections
 Bs = anuga.Transmissive_stage_zero_momentum_boundary(domain) # Strong reflections
 Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
@@ -100 +100 @@
     return -amplitude*sin((1./wave_length)*t*2*pi)
 
+def waveform2(t):
+    return 0.
+
 Bw2 = anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, waveform)
+Bf_in = swb2_boundary_conditions.Characteristic_boundary(domain, waveform)
 #Bw3 = swb2_boundary_conditions.Transmissive_momentum_nudge_stage_boundary(domain, waveform)
+Bf=swb2_boundary_conditions.Characteristic_boundary(domain,waveform2)
                     
 
 
 # Associate boundary tags with boundary objects
-domain.set_boundary({'left': Bw2, 'right': Bt, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Bw2, 'right': Bf, 'top': Br, 'bottom': Br})