Changeset 9292

Timestamp:

Aug 12, 2014, 12:10:51 PM (11 years ago)

Author:

davies

Message:

Adding local extrapolation and flux updating

Location:

trunk

Files:

• anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py (modified) (2 diffs)
• anuga_core/source/anuga/advection/advection.py (modified) (1 diff)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (8 diffs)
• anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c (modified) (35 diffs)
• anuga_core/source/anuga/shallow_water/test_local_extrapolation_and_flux_updating.py (added)
• anuga_work/development/gareth/tests/dam_break/dam_break.py (modified) (3 diffs)
• anuga_work/development/gareth/tests/dam_break/plotme.py (modified) (1 diff)
• anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py (modified) (3 diffs)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

@@ -1622 +1622 @@
 
             self.number_of_steps += 1
+
             if self._order_ == 1:
                 self.number_of_first_order_steps += 1
@@ -1678 +1679 @@
         # Update timestep to fit yieldstep and finaltime
         self.update_timestep(yieldstep, finaltime)
+
+        if self.max_flux_update_frequency is not 1:
+            # Update flux_update_frequency using the new timestep
+            self.compute_flux_update_frequency()
 
         # Update conserved quantities

trunk/anuga_core/source/anuga/advection/advection.py

@@ -80 +80 @@
         
         self.smooth = True
+        self.max_flux_update_frequency=1
 
     def check_integrity(self):

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

@@ -291 +291 @@
         self.boundary_flux_integral=boundary_flux_integral_operator(self)
         # Make an integer counting how many times we call compute_fluxes_central -- so we know which substep we are on
-        self.call=1 
+        #self.call=1 
 
         # List to store the volumes we computed before
@@ -297 +297 @@
         
         # Work arrays [avoid allocate statements in compute_fluxes or extrapolate_second_order]
-        self.edge_flux_work=num.zeros(len(self.edge_coordinates[:,0])*3)
-        self.pressuregrad_work=num.zeros(len(self.edge_coordinates[:,0]))
-        self.x_centroid_work=num.zeros(len(self.edge_coordinates[:,0])/3)
+        self.edge_flux_work=num.zeros(len(self.edge_coordinates[:,0])*3) # Advective fluxes
+        self.pressuregrad_work=num.zeros(len(self.edge_coordinates[:,0])) # Gravity related terms
+        self.x_centroid_work=num.zeros(len(self.edge_coordinates[:,0])/3) 
         self.y_centroid_work=num.zeros(len(self.edge_coordinates[:,0])/3)
+
+        ############################################################################
+        ## Local-timestepping information
+        #
+        # Fluxes can be updated every 1, 2, 4, 8, .. max_flux_update_frequency timesteps
+        # The global timestep is not allowed to increase except when
+        # number_of_timesteps%max_flux_update_frequency==0 
+        self.max_flux_update_frequency=2**0 # Must be a power of 2. 
+        
+        # flux_update_frequency. The edge flux terms are re-computed only when
+        #    number_of_timesteps%flux_update_frequency[myEdge]==0
+        self.flux_update_frequency=num.zeros(len(self.edge_coordinates[:,0])).astype(int)+1
+        # Flag: should we update the flux on the next compute fluxes call?
+        self.update_next_flux=num.zeros(len(self.edge_coordinates[:,0])).astype(int)+1
+        # Flag: should we update the extrapolation on the next extrapolation call?
+        # (Only do this if one or more of the fluxes on that triangle will be computed on
+        # the next timestep, assuming only the flux computation uses edge/vertex values)
+        self.update_extrapolation=num.zeros(len(self.edge_coordinates[:,0])/3).astype(int)+1
+
+        # edge_timestep [wavespeed/radius] -- not updated every timestep
+        self.edge_timestep=num.zeros(len(self.edge_coordinates[:,0]))+1.0e+100
+      
+        # Do we allow the timestep to increase (not every time if local
+        # extrapolation/flux updating is used) 
+        self.allow_timestep_increase=num.zeros(1).astype(int)+1
 
     def _set_defaults(self):
@@ -841 +866 @@
 
 
-
+    def set_local_extrapolation_and_flux_updating(self,nlevels=8):
+        """
+            Use local flux and extrapolation updating
+
+            nlevels == number of flux_update_frequency levels > 1
+
+                   For example, to allow flux updating every 1,2,4,8
+                   timesteps, do:
+
+                    domain.set_local_extrapolation_and_flux_updating(nlevels=3)
+
+                   (since 2**3==8)
+        """
+
+        self.max_flux_update_frequency=2**nlevels
+
+        if(self.max_flux_update_frequency is not 1):
+            if self.timestepping_method is not 'euler':
+                raise Exception, 'Local extrapolation and flux updating only supported with euler timestepping'
+            if self.compute_fluxes_method is not 'DE':
+                raise Exception, 'Local extrapolation and flux updating only supported for discontinuous flow algorithms'
 
 
@@ -1668 +1713 @@
             height = self.quantities['height']
 
-            timestep = float(sys.maxint)
-
-            # Keep track of number of calls to compute_fluxes_ext
-            # Note that in ANUGA, all sub-steps within an rk2/rk3 timestep have the same timestep
-            self.call+=1
+            timestep = self.evolve_max_timestep 
+
             flux_timestep = compute_fluxes_ext(timestep,
                                                self.epsilon,
@@ -1706 +1748 @@
                                                Bed.centroid_values, 
                                                height.centroid_values,
-                                               #Bed.vertex_values,
                                                self.edge_flux_type,
                                                self.riverwallData.riverwall_elevation,
@@ -1713 +1754 @@
                                                self.riverwallData.hydraulic_properties,
                                                self.edge_flux_work,
-                                               self.pressuregrad_work)
+                                               self.pressuregrad_work,
+                                               self.edge_timestep,
+                                               self.update_next_flux,
+                                               self.allow_timestep_increase)
 
             self.flux_timestep = flux_timestep
@@ -1821 +1865 @@
                       int(self.extrapolate_velocity_second_order),
                       self.x_centroid_work,
-                      self.y_centroid_work)
+                      self.y_centroid_work,
+                      self.update_extrapolation)
         else:
             # Code for original method
@@ -2580 +2625 @@
         return message        
           
+    def compute_flux_update_frequency(self):
+        """
+            Update the 'flux_update_frequency' and 'update_extrapolate' variables 
+            Used to control updating of fluxes / extrapolation for 'local-time-stepping'
+        """
+        from swDE1_domain_ext import compute_flux_update_frequency \
+                                  as compute_flux_update_frequency_ext
+
+        compute_flux_update_frequency_ext(
+                self.timestep,
+                self.neighbours,
+                self.neighbour_edges,
+                self.already_computed_flux,
+                self.edge_timestep,
+                self.flux_update_frequency,
+                self.update_extrapolation,
+                self.max_flux_update_frequency,
+                self.update_next_flux,
+                self.allow_timestep_increase)
+        
     def report_water_volume_statistics(self, verbose=True, returnStats=False):
         """

trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c

@@ -28 +28 @@
 
 const double pi = 3.14159265358979;
+
+// Trick to compute n modulo 2 when d is a power of 2
+unsigned int Mod_of_power_2(unsigned int n, unsigned int d)
+{
+  return ( n & (d-1) );
+} 
+
 
 // Computational function for rotation
@@ -111 +118 @@
                            double *edgeflux, double *max_speed,
                            double *pressure_flux, double hc,
-                           double hc_n, double speed_max_last)
+                           double hc_n)
 {
 
@@ -285 +292 @@
                            double *edgeflux, double *max_speed, 
                            double *pressure_flux, double hc, 
-                           double hc_n, double speed_max_last) 
+                           double hc_n) 
 {
 
@@ -441 +448 @@
   return 0;
 }
+
+////////////////////////////////////////////////////////////////
+
+int _compute_flux_update_frequency(int number_of_elements, 
+                                  long *neighbours, long *neighbour_edges,
+                                  long *already_computed_flux,
+                                  int max_flux_update_frequency, double *edge_timestep,
+                                  double timestep, double notSoFast,
+                                  long *flux_update_frequency,
+                                  long *update_extrapolation, 
+                                  long *update_next_flux,
+                                  long *allow_timestep_increase){
+    // Compute the 'flux_update_frequency' for each edge.
+    //
+    // This determines how regularly we need
+    // to update the flux computation (not every timestep)
+    //
+    // Allowed values are 1,2,4,8,... max_flux_update_frequency.
+    //
+    // For example, an edge with flux_update_frequency = 4 would
+    // only have the flux updated every 4 timesteps
+    //
+    //
+    // Local variables
+    int k, i, k3, ki, m, n, nm, ii, j, ii2;
+    long fuf;
+    static int cyclic_number_of_steps=-1;
+
+    // QUICK EXIT
+    if(max_flux_update_frequency==1){
+        return 0;
+    }
+    
+    // Count the steps
+    cyclic_number_of_steps++;
+    if(cyclic_number_of_steps==max_flux_update_frequency){
+        // The flux was just updated in every cell
+        cyclic_number_of_steps=0;
+    }
+
+
+    // PART 1: ONLY OCCURS FOLLOWING FLUX UPDATE
+
+    for ( k = 0; k < number_of_elements; k++){
+        for ( i = 0; i < 3; i++){
+            ki = k*3 + i;
+            if((Mod_of_power_2(cyclic_number_of_steps, flux_update_frequency[ki])==0)){
+                // The flux was just updated, along with the edge_timestep
+                // So we better recompute the flux_update_frequency
+                n=neighbours[ki];
+                if(n>=0){
+                    m = neighbour_edges[ki];
+                    nm = n * 3 + m; // Linear index (triangle n, edge m)
+                }
+
+                // Check if we have already done this edge
+                // (Multiply already_computed_flux by -1 on the first update,
+                // and again on the 2nd)
+                if(already_computed_flux[ki] > 0 ){
+                    // We have not fixed this flux value yet
+                    already_computed_flux[ki] *=-1;
+                    if(n>=0){
+                        already_computed_flux[nm] *=-1;
+                    }
+                }else{
+                    // We have fixed this flux value already 
+                    already_computed_flux[ki] *=-1;
+                    if(n>=0){
+                        already_computed_flux[nm] *=-1;
+                    }
+                    continue;
+                }
+
+                // Basically int( edge_ki_timestep/timestep ) with upper limit + tweaks
+                // notSoFast is ideally = 1.0, but in practice values < 1.0 can enhance stability
+                // NOTE: edge_timestep[ki]/timestep can be very large [so int overflows].
+                //       Do not pull the (int) inside the min term
+                fuf = (int)min((edge_timestep[ki]/timestep)*notSoFast,max_flux_update_frequency*1.);
+                // Account for neighbour
+                if(n>=0){
+                    fuf = min( (int)min(edge_timestep[nm]/timestep*notSoFast, max_flux_update_frequency*1.), fuf);
+                }
+
+                // Deal with notSoFast<1.0
+                if(fuf<1){
+                    fuf=1;
+                }
+                // Deal with large fuf 
+                if(fuf>max_flux_update_frequency){
+                    fuf = max_flux_update_frequency;
+                }
+                //// Deal with intermediate cases
+                ii=2;
+                while(ii<max_flux_update_frequency){
+                    // Set it to 1,2,4, 8, ...
+                    ii2=2*ii;
+                    if((fuf>ii) && (fuf<ii2)){
+                        fuf = ii;
+                        continue;
+                    }
+                    ii=ii2;
+                }
+
+                // Set the values
+                flux_update_frequency[ki]=fuf;
+                if(n>=0){
+                    flux_update_frequency[nm]= fuf;
+                }
+                
+            }
+        }
+    } 
+
+    //// PART 2 -- occcurs every timestep
+
+    // At this point, both edges have the same flux_update_frequency.
+    // Next, ensure that flux_update_frequency varies within a constant over each triangle
+    // Experiences suggests this is numerically important
+    // (But, it can result in the same edge having different flux_update_freq)
+    for( k=0; k<number_of_elements; k++){
+        k3=3*k;
+        ii = 1*min(flux_update_frequency[k3],
+                 min(flux_update_frequency[k3+1],
+                     flux_update_frequency[k3+2]));
+        
+        flux_update_frequency[k3]=min(ii, flux_update_frequency[k3]);
+        flux_update_frequency[k3+1]=min(ii, flux_update_frequency[k3+1]);
+        flux_update_frequency[k3+2]=min(ii,flux_update_frequency[k3+2]);
+            
+    }
+            
+    // Now enforce the same flux_update_frequency on each edge
+    // (Could have been broken above when we limited the variation on each triangle)
+    // This seems to have nice behaviour. Notice how an edge
+    // with a large flux_update_frequency, near an edge with a small flux_update_frequency,
+    // will have its flux_update_frequency updated after a few timesteps (i.e. before max_flux_update_frequency timesteps)
+    // OTOH, could this cause oscillations in flux_update_frequency?
+    for( k = 0; k < number_of_elements; k++){
+        update_extrapolation[k]=0;
+        for( i = 0; i< 3; i++){
+            ki=3*k+i;
+            // Account for neighbour
+            n=neighbours[ki];
+            if(n>=0){
+                m = neighbour_edges[ki];
+                nm = n * 3 + m; // Linear index (triangle n, edge m)
+                flux_update_frequency[ki]=min(flux_update_frequency[ki], flux_update_frequency[nm]);
+            }
+            // Do we need to update the extrapolation?
+            // (We do if the next flux computation will actually compute a flux!)
+            if(Mod_of_power_2((cyclic_number_of_steps+1),flux_update_frequency[ki])==0){
+                update_next_flux[ki]=1;
+                update_extrapolation[k]=1;
+            }else{
+                update_next_flux[ki]=0;
+            }
+        }
+    }
+
+    // Check whether the timestep can be increased in the next compute_fluxes call
+    if(cyclic_number_of_steps+1==max_flux_update_frequency){
+        // All fluxes will be updated on the next timestep
+        // We also allow the timestep to increase then
+        allow_timestep_increase[0]=1;
+    }else{
+        allow_timestep_increase[0]=0;
+    }
+
+    return 0;
+}
+
 
 double adjust_edgeflux_with_weir(double* edgeflux,
@@ -574 +752 @@
         double* riverwall_hydraulic_properties,
         double* edge_flux_work,
-        double* pressuregrad_work) {
+        double* pressuregrad_work,
+        double* edge_timestep,
+        long* update_next_flux,
+        long* allow_timestep_increase) {
     // Local variables
     double max_speed, length, inv_area, zl, zr;
@@ -582 +763 @@
     //
     int k, i, m, n,j, ii;
-    int ki, nm = 0, ki2,ki3, nm3; // Index shorthands
+    int ki,k3, nm = 0, ki2,ki3, nm3; // Index shorthands
     // Workspace (making them static actually made function slightly slower (Ole))
     double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
     double stage_edges[3];//Work array
-    double bedslope_work, local_timestep;
+    double bedslope_work;
+    static double local_timestep;
     int neighbours_wet[3];//Work array
-    int RiverWall_count, substep_count;
+    long RiverWall_count, substep_count;
     double hle, hre, zc, zc_n, Qfactor, s1, s2, h1, h2; 
-    double stage_edge_lim, outgoing_mass_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2;
+    double stage_edge_lim, outgoing_mass_edges, pressure_flux, hc, hc_n, tmp, tmp2;
     double h_left_tmp, h_right_tmp;
     static long call = 1; // Static local variable flagging already computed flux
-    double speed_max_last, vol, smooth, local_speed, weir_height;
+    double speed_max_last, vol, weir_height;
 
     call++; // Flag 'id' of flux calculation for this timestep
@@ -602 +784 @@
     memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double));
     memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
-    memset((char*) edge_flux_work, 0, 9*number_of_elements * sizeof (double));
-    memset((char*) pressuregrad_work, 0, 3*number_of_elements * sizeof (double));
-
-    local_timestep=timestep;
+
+
+    // Counter for riverwall edges
     RiverWall_count=0;
+    // Which substep of the timestepping method are we on?
     substep_count=(call-2)%timestep_fluxcalls;
-
+    
+    // Fluxes are not updated every timestep,
+    // but all fluxes ARE updated when the following condition holds
+    if(allow_timestep_increase[0]==1){
+        // We can only increase the timestep if all fluxes are allowed to be updated
+        // If this is not done the timestep can't increase (since local_timestep is static)
+        local_timestep=1.0e+100;
+    }
 
     // For all triangles
     for (k = 0; k < number_of_elements; k++) {
         speed_max_last=0.0;
+
         // Loop through neighbours and compute edge flux for each
         for (i = 0; i < 3; i++) {
@@ -619 +809 @@
             ki3 = 3*ki; 
 
-            if (already_computed_flux[ki] == call) {
+            if ((already_computed_flux[ki] == call) || (update_next_flux[ki]!=1)) {
                 // We've already computed the flux across this edge
                 // Check if it is a riverwall
@@ -693 +883 @@
                     normals[ki2], normals[ki2 + 1],
                     epsilon, z_half, limiting_threshold, g,
-                    edgeflux, &max_speed, &pressure_flux, hc, hc_n, 
-                    speed_max_last);
+                    edgeflux, &max_speed, &pressure_flux, hc, hc_n);
 
             // Force weir discharge to match weir theory
@@ -701 +890 @@
                 // If the weir height is zero, avoid the weir computation entirely
                 if(weir_height>0.){
+                    ////////////////////////////////////////////////////////////////////////////////////
+                    // Use first-order h's for weir -- as the 'upstream/downstream' heads are
+                    //  measured away from the weir itself
+                    h_left_tmp= max(stage_centroid_values[k]-z_half,0.);
+                    if(n>=0){
+                        h_right_tmp= max(stage_centroid_values[n]-z_half,0.);
+                    }else{
+                        h_right_tmp= max(hc_n+zr-z_half,0.);
+                    }
+
+                    //////////////////////////////////////////////////////////////////////////////////
                     // Get Qfactor index - multiply the idealised weir discharge by this constant factor
                     ii=riverwall_rowIndex[RiverWall_count-1]*ncol_riverwall_hydraulic_properties;
@@ -717 +917 @@
                     h2=riverwall_hydraulic_properties[ii];
                     
-                    // Use first-order h's for weir -- as the 'upstream/downstream' heads are
-                    //  measured away from the weir itself
-                    h_left_tmp= max(stage_centroid_values[k]-z_half,0.);
-                    if(n>=0){
-                        h_right_tmp= max(stage_centroid_values[n]-z_half,0.);
-                    }else{
-                        h_right_tmp= max(hc_n+zr-z_half,0.);
-                    }
-
+                    // Weir flux adjustment 
                     adjust_edgeflux_with_weir(edgeflux, h_left_tmp, h_right_tmp, g, 
                                               weir_height, Qfactor, 
@@ -782 +974 @@
 
             // Update timestep based on edge i and possibly neighbour n
-            // NOTE: We should only change the timestep between rk2 or rk3
-            // steps, NOT within them (since a constant timestep is used within
-            // each rk2/rk3 sub-step)
-            if ((tri_full_flag[k] == 1) & (substep_count==0)) {
-
-                speed_max_last=max(speed_max_last, max_speed);
-
-                if (max_speed > epsilon) {
-                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
-
-                    // CFL for triangle k
-                    local_timestep = min(local_timestep, radii[k] / max_speed);
-
-                    if (n >= 0) {
-                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
-                        local_timestep = min(local_timestep, radii[n] / max_speed);
+            // NOTE: We should only change the timestep on the 'first substep'
+            //  of the timestepping method [substep_count==0]
+            if(substep_count==0){
+
+                // Compute the 'edge-timesteps' (useful for setting flux_update_frequency)
+                edge_timestep[ki] = radii[k] / max(max_speed, epsilon);
+                if (n >= 0) {
+                    edge_timestep[nm] = radii[n] / max(max_speed, epsilon);
+                }
+
+                // Update the timestep
+                if ((tri_full_flag[k] == 1)) {
+
+                    speed_max_last=max(speed_max_last, max_speed);
+
+                    if (max_speed > epsilon) {
+                        // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
+
+                        // CFL for triangle k
+                        local_timestep = min(local_timestep, edge_timestep[ki]);
+
+                        if (n >= 0) {
+                            // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
+                            local_timestep = min(local_timestep, edge_timestep[nm]);
+                        }
                     }
-
                 }
             }
-        
+
         } // End edge i (and neighbour n)
         // Keep track of maximal speeds
@@ -809 +1009 @@
 
     } // End triangle k
- 
+
     //// Limit edgefluxes, for mass conservation near wet/dry cells
     //// This doesn't seem to be needed anymore
@@ -867 +1067 @@
     for(k=0; k<number_of_elements; k++){
         hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.);
-        stage_explicit_update[k]=0.;
-        xmom_explicit_update[k]=0.;
-        ymom_explicit_update[k]=0.;
 
         for(i=0;i<3;i++){
+            // FIXME: Make use of neighbours to efficiently set things
             ki=3*k+i;   
             ki2=ki*2;
@@ -906 +1104 @@
             //    ymom_explicit_update[k] *= 0.;
             //}
+            
 
         } // end edge i
@@ -920 +1119 @@
     // Ensure we only update the timestep on the first call within each rk2/rk3 step
     if(substep_count==0) timestep=local_timestep; 
-            
+         
     return timestep;
 }
@@ -1075 +1274 @@
                                  int extrapolate_velocity_second_order,
                                  double* x_centroid_work,
-                                 double* y_centroid_work) {
+                                 double* y_centroid_work,
+                                 long* update_extrapolation) {
                   
   // Local variables
@@ -1122 +1322 @@
   // some situations
   for (k=0; k<number_of_elements;k++){
-
+      
       k3=k*3;
       k0 = surrogate_neighbours[k3];
@@ -1145 +1345 @@
   for (k = 0; k < number_of_elements; k++) 
   {
+
+    // Don't update the extrapolation if the flux will not be computed on the
+    // next timestep
+    if(update_extrapolation[k]==0){
+       continue;
+    }
+
 
     // Useful indices
@@ -1223 +1430 @@
       k2 = surrogate_neighbours[k3 + 2];
 
+
       // Get the auxiliary triangle's vertex coordinates 
       // (really the centroids of neighbouring triangles)
@@ -1264 +1472 @@
           height_edge_values[k3+2] = dk;
           
-          //stage_edge_values[k3]   = elevation_centroid_values[k]+dk;
-          //stage_edge_values[k3+1] = elevation_centroid_values[k]+dk;
-          //stage_edge_values[k3+2] = elevation_centroid_values[k]+dk;
-          //stage_edge_values[k3] = elevation_edge_values[k3]+dk;
-          //stage_edge_values[k3+1] = elevation_edge_values[k3+1]+dk;
-          //stage_edge_values[k3+2] = elevation_edge_values[k3+2]+dk;
-
-          //xmom_centroid_values[k]=0.;
-          //ymom_centroid_values[k]=0.;
-          //x_centroid_work[k] = 0.;
-          //y_centroid_work[k] = 0.;
-
           xmom_edge_values[k3]    = xmom_centroid_values[k];
           xmom_edge_values[k3+1]  = xmom_centroid_values[k];
@@ -1711 +1907 @@
   // Compute vertex values of quantities
   for (k=0; k<number_of_elements; k++){
+      if(extrapolate_velocity_second_order==1){
+          //Convert velocity back to momenta at centroids
+          xmom_centroid_values[k] = x_centroid_work[k];
+          ymom_centroid_values[k] = y_centroid_work[k];
+      }
+     
+      // Don't proceed if we didn't update the edge/vertex values
+      if(update_extrapolation[k]==0){
+         continue;
+      }
+
       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] ; 
@@ -1725 +1932 @@
       // If needed, convert from velocity to momenta
       if(extrapolate_velocity_second_order==1){
-          //Convert velocity back to momenta at centroids
-          xmom_centroid_values[k] = x_centroid_work[k];
-          ymom_centroid_values[k] = y_centroid_work[k];
-      
           // Re-compute momenta at edges
           for (i=0; i<3; i++){
@@ -1847 +2050 @@
     *riverwall_hydraulic_properties,
     *edge_flux_work,
-    *pressuregrad_work;
+    *pressuregrad_work,
+    *edge_timestep,
+    *update_next_flux,
+    *allow_timestep_increase;
     
   double timestep, epsilon, H0, g;
   int optimise_dry_cells, timestep_fluxcalls, ncol_riverwall_hydraulic_properties;
+  
     
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOOiiOOOOOOOOiOOO",
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOOiiOOOOOOOOiOOOOOO",
             &timestep,
             &epsilon,
@@ -1892 +2099 @@
             &riverwall_hydraulic_properties,
             &edge_flux_work,
-            &pressuregrad_work
+            &pressuregrad_work,
+            &edge_timestep,
+            &update_next_flux,
+            &allow_timestep_increase
             )) {
     report_python_error(AT, "could not parse input arguments");
@@ -1932 +2142 @@
   CHECK_C_CONTIG(edge_flux_work);
   CHECK_C_CONTIG(pressuregrad_work);
+  CHECK_C_CONTIG(edge_timestep);
+  CHECK_C_CONTIG(update_next_flux);
+  CHECK_C_CONTIG(allow_timestep_increase);
 
   int number_of_elements = stage_edge_values -> dimensions[0];
@@ -1978 +2191 @@
                      (double*) riverwall_hydraulic_properties ->data, 
                      (double*) edge_flux_work-> data,
-                     (double*) pressuregrad_work->data);
+                     (double*) pressuregrad_work->data,
+                     (double*) edge_timestep->data,
+                     (long*) update_next_flux->data,
+                     (long*) allow_timestep_increase->data
+                     );
   // Return updated flux timestep
   return Py_BuildValue("d", timestep);
@@ -2029 +2246 @@
              &pressure_flux,
                          ((double*) normal -> data)[0],
-                         ((double*) normal -> data)[1], 
                          ((double*) normal -> data)[1]
              );
@@ -2107 +2323 @@
 }
 
+PyObject *compute_flux_update_frequency(PyObject *self, PyObject *args) {
+  /*Compute how often we should update fluxes and extrapolations (not every timestep)
+
+    python_call
+        compute_flux_update_frequency_ext(
+                self.timestep,
+                self.neighbours,
+                self.neighbour_edges,
+                self.already_computed_flux,
+                self.edge_timestep,
+                self.flux_update_frequency,
+                self.update_extrapolation,
+                self.max_flux_update_frequency,
+                self.update_next_flux)
+
+
+
+  */
+
+
+  PyArrayObject *neighbours, *neighbour_edges,
+    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
+    *edge_timestep,
+    *flux_update_frequency,
+    *update_extrapolation,
+    *update_next_flux,
+    *allow_timestep_increase;
+    
+  double timestep;
+  int max_flux_update_frequency;
+  
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "dOOOOOOiOO",
+            &timestep,
+            &neighbours,
+            &neighbour_edges,
+            &already_computed_flux,
+            &edge_timestep,
+            &flux_update_frequency,
+            &update_extrapolation,
+            &max_flux_update_frequency,
+            &update_next_flux,
+            &allow_timestep_increase
+            )) {
+    report_python_error(AT, "could not parse input arguments");
+    return NULL;
+  }
+
+  // check that numpy array objects arrays are C contiguous memory
+  CHECK_C_CONTIG(neighbours);
+  CHECK_C_CONTIG(neighbour_edges);
+  CHECK_C_CONTIG(already_computed_flux);
+  CHECK_C_CONTIG(edge_timestep);
+  CHECK_C_CONTIG(flux_update_frequency);
+  CHECK_C_CONTIG(update_extrapolation);
+  CHECK_C_CONTIG(update_next_flux);
+  CHECK_C_CONTIG(allow_timestep_increase);
+
+  int number_of_elements = update_extrapolation -> dimensions[0];
+  // flux_update_frequency determined by rounding edge_timestep/timestep*notSoFast
+  // So notSoFast<1 might increase the number of flux calls a cell has to do, but
+  // can be useful for increasing numerical stability
+  double notSoFast=1.00;
+
+  // Call underlying flux computation routine and update 
+  // the explicit update arrays 
+  _compute_flux_update_frequency(number_of_elements,
+                                 (long*) neighbours->data,
+                                 (long*) neighbour_edges->data,
+                                 (long*) already_computed_flux->data,
+                                 max_flux_update_frequency,
+                                 (double*) edge_timestep->data,
+                                 timestep,
+                                 notSoFast,
+                                 (long*) flux_update_frequency->data,
+                                 (long*) update_extrapolation->data,
+                                 (long*) update_next_flux->data,
+                                 (long*) allow_timestep_increase->data
+                                 );
+  // Return 
+  return Py_BuildValue("");
+}
+
 
 PyObject *extrapolate_second_order_edge_sw(PyObject *self, PyObject *args) {
@@ -2151 +2451 @@
     *height_vertex_values,
     *x_centroid_work,
-    *y_centroid_work;
+    *y_centroid_work,
+    *update_extrapolation;
   
   PyObject *domain;
@@ -2161 +2462 @@
   
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOOOOOiiOO",
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOOOOOiiOOO",
             &domain,
             &surrogate_neighbours,
@@ -2186 +2487 @@
             &extrapolate_velocity_second_order,
             &x_centroid_work,
-            &y_centroid_work)) {         
+            &y_centroid_work,
+            &update_extrapolation)) {         
 
     report_python_error(AT, "could not parse input arguments");
@@ -2215 +2517 @@
   CHECK_C_CONTIG(x_centroid_work);
   CHECK_C_CONTIG(y_centroid_work);
+  CHECK_C_CONTIG(update_extrapolation);
   
   // Get the safety factor beta_w, set in the config.py file. 
@@ -2267 +2570 @@
                    extrapolate_velocity_second_order,
                    (double*) x_centroid_work -> data,
-                   (double*) y_centroid_work -> data);
+                   (double*) y_centroid_work -> data,
+                   (long*) update_extrapolation -> data);
 
   if (e == -1) {
@@ -2346 +2650 @@
   {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},  
   {"extrapolate_second_order_edge_sw", extrapolate_second_order_edge_sw, METH_VARARGS, "Print out"},
+  {"compute_flux_update_frequency", compute_flux_update_frequency, METH_VARARGS, "Print out"},
   {"protect",          protect, METH_VARARGS | METH_KEYWORDS, "Print out"},
   {NULL, NULL, 0, NULL}

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

@@ -14 +14 @@
 from numpy import zeros, float
 from time import localtime, strftime, gmtime
-from bal_and import *
+#from bal_and import *
 #from anuga_tsunami import *
 
@@ -42 +42 @@
 points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
 
-#domain = anuga.Domain(points, vertices, boundary) 
-domain = Domain(points, vertices, boundary) 
+domain = anuga.Domain(points, vertices, boundary) 
+#domain = Domain(points, vertices, boundary) 
 
 domain.set_name(output_file)                
@@ -56 +56 @@
 print domain.get_timestepping_method()
 
+domain.set_flow_algorithm('DE0')
+#domain.set_local_extrapolation_and_flux_updating(nlevels=8)
 #domain.use_edge_limiter = True
 #domain.use_edge_limiter = False

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

@@ -4 +4 @@
 """
 from anuga.utilities import plot_utils as util
-from bal_and import plot_utils as util2
 from matplotlib import pyplot as pyplot
 
-p_st = util.get_output('dam_break_20120305_124724/dam_break.sww')
+p_st = util.get_output('dam_break_20140811_130209/dam_break.sww')
 p2_st=util.get_centroids(p_st)
 
 
-p_dev = util2.get_output('dam_break_20131204_145613/dam_break.sww', 0.001)
-p2_dev=util2.get_centroids(p_dev, velocity_extrapolation=True)
+p_dev = util.get_output('dam_break_20140811_130049/dam_break.sww', 0.001)
+p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
 
 pyplot.clf()

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

@@ -20 +20 @@
 domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
 #domain.set_store_centroids(True)
-domain.set_flow_algorithm('DE1')
+domain.set_flow_algorithm('DE0')
 domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 
          'ymomentum': 2, 'elevation': 2})
@@ -42 +42 @@
 bumpiness=50. # Higher = shorter wavelength oscillations in topography
 tstep=0.2
-lasttime=90.
+lasttime=30.
 
 #domain.minimum_allowed_height=domain.minimum_allowed_height*scale_me # Seems needed to make the algorithms behave
@@ -110 +110 @@
     print 'Their difference is:', vol-flux_integral
     #print 'Volume less flux int', sum(dd_raw*domain.areas) - domain.boundary_flux_integral
-    
+    print 'Update_next_flux range:'
+    print domain.update_extrapolation.min(), domain.update_extrapolation.max()
+    print 'Flux_update_frequency range:'
+    print domain.flux_update_frequency.min(), domain.flux_update_frequency.max()
 
 vel_final_inds=(vv>1.0e-01).nonzero()