Changeset 8751

Timestamp:

Mar 11, 2013, 8:17:49 PM (12 years ago)

Author:

davies

Message:

Updates to gd in-development algorithms

Location:

trunk

Files:

• anuga_core/source/anuga_validation_tests/utilities/fabricate.py (modified) (1 prop)
• anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py (modified) (3 diffs)
• anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c (modified) (37 diffs)
• anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py (modified) (1 diff)
• anuga_work/development/gareth/tests/dam_break/dam_break.py (modified) (1 diff)
• anuga_work/development/gareth/tests/dam_break/plotme.py (modified) (1 diff)
• anuga_work/development/gareth/tests/okushiri/run_okushiri.py (modified) (1 diff)
• anuga_work/development/gareth/tests/parabolic/parabolaplot.py (modified) (1 diff)
• anuga_work/development/gareth/tests/runup/runup.py (modified) (1 diff)
• anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py (modified) (8 diffs)
• anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py (modified) (1 diff)
• anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py (modified) (2 diffs)
• anuga_work/development/gareth/tests/wave_runup (added)
• anuga_work/development/gareth/tests/wave_runup/Shoreline.csv (added)
• anuga_work/development/gareth/tests/wave_runup/initial_condition.txt (added)
• anuga_work/development/gareth/tests/wave_runup/plotme.py (added)
• anuga_work/development/gareth/tests/wave_runup/readme.txt (added)
• anuga_work/development/gareth/tests/wave_runup/runup.py (added)
• anuga_work/development/gareth/tests/wave_runup/t160.csv (added)
• anuga_work/development/gareth/tests/wave_runup/t175.csv (added)
• anuga_work/development/gareth/tests/wave_runup/t220.csv (added)
• anuga_work/development/gareth/tests/wave_runup/util.py (added)

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

@@ -63 +63 @@
         # Most of these override the options in config.py
         #------------------------------------------------
-        self.set_CFL(1.0)
+        self.set_CFL(1.00)
         self.set_use_kinematic_viscosity(False)
-        self.timestepping_method='rk2' 
+        self.timestepping_method='rk2' #'euler'#'rk2'#'euler'#'rk2' 
         # The following allows storage of the negative depths associated with this method
         self.minimum_storable_height=-99999999999.0 
@@ -235 +235 @@
 
         # Shortcuts
+        if(domain.timestepping_method!='rk2'):
+            err_mess='ERROR: Timestepping method is ' + domain.timestepping_method +'. '+\
+                     'You need to use rk2 timestepping with this solver, ' +\
+                     ' because there is presently a hack in compute fluxes central. \n'+\
+                     ' The HACK: The timestep will only ' +\
+                     'be recomputed on every 2nd call to compute_fluxes_central, to support'+\
+                     ' correct treatment of wetting and drying'
+
+            raise Exception, err_mess 
+
         Stage = domain.quantities['stage']
         Xmom = domain.quantities['xmomentum']
@@ -268 +278 @@
                                            int(domain.optimise_dry_cells),
                                            Stage.centroid_values,
+                                           Xmom.centroid_values,
+                                           Ymom.centroid_values,
                                            Bed.centroid_values, 
                                            Bed.vertex_values)

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

@@ -206 +206 @@
   denom = s_max - s_min;
   if (denom < epsilon) 
-  { // FIXME (Ole): Try using h0 here
+  { 
+    // Both wave speeds are very small
     memset(edgeflux, 0, 3*sizeof(double));
 
@@ -263 +264 @@
         int optimise_dry_cells, 
         double* stage_centroid_values,
+        double* xmom_centroid_values,
+        double* ymom_centroid_values,
         double* bed_centroid_values,
         double* bed_vertex_values) {
@@ -273 +276 @@
     // See ANUGA manual under flux limiting
     int k, i, m, n,j;
-    int ki, nm = 0, ki2; // Index shorthands
+    int ki, nm = 0, ki2,ki3, nm3; // Index shorthands
 
     // Workspace (making them static actually made function slightly slower (Ole))
@@ -281 +284 @@
     int neighbours_wet[3];//Work array
     int useint;
-    double stage_edge_lim, scale_factor_shallow, bedtop, bedbot, pressure_flux, hc, hc_n, tmp;
-    //double max_bed_edgevalue[number_of_elements]; // FIXME: Set memory explicitly
+    double stage_edge_lim, outgoing_flux_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2;
     static long call = 1; // Static local variable flagging already computed flux
 
-    double *max_bed_edgevalue, *min_bed_edgevalue;
+    double *max_bed_edgevalue, *min_bed_edgevalue, *edgeflux_store, *pressuregrad_store;
 
     max_bed_edgevalue = malloc(number_of_elements*sizeof(double));
-    //min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
-    // Start computation
+    edgeflux_store = malloc(number_of_elements*9*sizeof(double));
+    pressuregrad_store = malloc(number_of_elements*3*sizeof(double));
+
     call++; // Flag 'id' of flux calculation for this timestep
 
@@ -298 +301 @@
     memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
 
-
     // Compute minimum 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]));
-        //max_bed_edgevalue[k]= 0.5*(max_bed_edgevalue[j]+bed_centroid_values[k]);
-        //min_bed_edgevalue[k] = min(bed_edge_values[3*k], 
-        //                           min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
-    
     }
 
@@ -315 +313 @@
         for (i = 0; i < 3; i++) {
             ki = k * 3 + i; // Linear index to edge i of triangle k
+            ki2 = 2 * ki; //k*6 + i*2
+            ki3 = 3*ki; 
 
             if (already_computed_flux[ki] == call) {
@@ -328 +328 @@
             hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.0);
 
-            //if(ql[0]<=zl){
-            //    ql[1]=0.;
-            //    ql[2]=0.;
-            //}
             // Get right hand side values either from neighbouring triangle
             // or from boundary array (Quantities at neighbour on nearest face).
@@ -350 +346 @@
                 m = neighbour_edges[ki];
                 nm = n * 3 + m; // Linear index (triangle n, edge m)
+                nm3 = nm*3;
 
                 qr[0] = stage_edge_values[nm];
@@ -357 +354 @@
             }
             
-
             if (fabs(zl-zr)>1.0e-10) {
                 report_python_error(AT,"Discontinuous Elevation");
                 return 0.0;
             }
-            
+           
+            // GD HACK -- although I think this is a common technique 
             //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
             // unless the local centroid value is smaller 
             if(n>=0){
-                //if(hc==0.0){
                 if((hc<=H0)&&(hc_n>H0)){
                     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]=max(min(qr[0],stage_centroid_values[k]),zl);
                 }
-                //if(hc_n==0.0){
                 if((hc_n<=H0)&&(hc>H0)){
                     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]; 
                 }
-
-                //if((hc_n<=H0)&&(hc<=H0)){
-                //    ql[0] = zl;
-                //    qr[0] = zr;
-                //}
 
             }else{
                 // Treat the boundary case
-                //if((hc==0.0)){
-                //    ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
-                //    //ql[0]=max(min(qr[0],ql[0]),zl);
-                //}
+                // ??
             }
-            //if(k==229 | n==229){
-            //    printf(" ---------------\n ");
-            //    printf("k, ql, qr, zl, zr==%d, %20.20e, %20.20e, %20.20e, %20.20e \n", k, ql[0], qr[0], zl, zr);
-
-            //}
-            
-            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
-            ki2 = 2 * ki; //k*6 + i*2
 
             // Edge flux computation (triangle k, edge i)
@@ -406 +381 @@
                     edgeflux, &max_speed, &pressure_flux);
 
-            // Prevent outflow from 'seriously' dry cells
-            // Idea: The cell will not go dry if:
-            // mass_flux = edgeflux[0]*(dt*edgelength) <= (1/3)Area_triangle*hc
-            if((stage_centroid_values[k]<=max_bed_edgevalue[k])|
-               (ql[0]<=zl)){
-                if(edgeflux[0]>0.0){
-                    tmp=min((1.0/3.0)*areas[k]*hc/(edgelengths[ki]*max(timestep,epsilon)), 1.0); // 28/7 -- Right answer for channel flow problem.
-                    tmp = min(min(edgeflux[0], tmp)/edgeflux[0], 1.0);
-                    if((tmp<0.)|tmp>1.0){
-                        printf("Invalid tmp value\n");
-                    }
-                    edgeflux[0]*=tmp;
-                    edgeflux[1]*=tmp;
-                    edgeflux[2]*=tmp;
-                }
-            }
-            if(n>=0){
-                if((stage_centroid_values[n]<=max_bed_edgevalue[n])|
-                   (qr[0]<=zr)){
-                    if(edgeflux[0]<0.0){
-                        tmp=min((1.0/3.0)*areas[n]*hc_n/(edgelengths[ki]*max(timestep,epsilon)), 1.0); // 28/7 -- Right answer for channel flow problem.
-                        tmp = min( max(edgeflux[0], -tmp)/edgeflux[0], 1.0);
-                        if((tmp<0.)|tmp>1.0){
-                            printf("Invalid tmp value\n");
-                        }
-                        edgeflux[0]*=tmp;
-                        edgeflux[1]*=tmp;
-                        edgeflux[2]*=tmp;
-                    }
-                }
-            } 
            
-            //if(k==229|n==229){
-            //    printf("++edgeflux=, %20.20e, %20.20e, %20.20e \n", edgeflux[0],edgeflux[1], edgeflux[2]);
-            //}
-
             // Multiply edgeflux by edgelength
             length = edgelengths[ki];
@@ -448 +388 @@
             edgeflux[2] *= length;
 
+            edgeflux_store[ki3 + 0 ] = -edgeflux[0];
+            edgeflux_store[ki3 + 1 ] = -edgeflux[1];
+            edgeflux_store[ki3 + 2 ] = -edgeflux[2];
 
             // Update triangle k with flux from edge i
-            stage_explicit_update[k] -= edgeflux[0];
-            xmom_explicit_update[k] -= edgeflux[1];
-            ymom_explicit_update[k] -= edgeflux[2];
-
-            // Compute bed slope term
-            //if(hc>-9999.0){
-                //Bedslope approx 1:
-            if(hc>H0){
-                // Add the pressure gradient flux term back in
-                // Do this to avoid the above mass-conservative limiting
-                // Because limiting the pressure term would destroy balancing with the bedslope term
-                //edgeflux[1]+=normals[ki2]*pressure_flux;
-                //edgeflux[2]+=normals[ki2+1]*pressure_flux;
-
-                //bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
-                bedslope_work = length*g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux*length;
-                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
-                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
-            }else{
-                //bedslope_work = 0.;//-pressure_flux*length ;//g*length*( 0.*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
-                xmom_explicit_update[k] *= 0.;
-                ymom_explicit_update[k] *= 0.;
-            }
+
+            // bedslope_work contains all gravity related terms -- weighting of
+            // conservative and non-conservative versions.
+            bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux);
+            //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 
+            pressuregrad_store[ki] = bedslope_work;
+
             
-            //if(k==229|n==229){
-            //    printf("--k-pres, xmom, bedwrk, %20.20e, %20.20e, %20.20e \n", pressure_flux*length, xmom_explicit_update[k], bedslope_work);
-            ////    printf(" Edgesum: %20.20e\n", edgelengths[3*k]*normals[6*k]*g*hc*stage_edge_values[3*k]
-            ////                   + edgelengths[3*k+1]*normals[6*k+2]*g*hc*stage_edge_values[3*k+1]
-            ////                   + edgelengths[3*k+2]*normals[6*k+4]*g*hc*stage_edge_values[3*k+2]);
-            //}
-                //
-                // Bedslope approx 2
-                //stage_edge_lim = ql[0]; // Limit this to be between a constant stage and constant depth extrapolation
-                //if(stage_edge_lim > max(stage_centroid_values[k], zl +hc)){
-                //    stage_edge_lim = max(stage_centroid_values[k], zl+hc);
-                //}
-                //if(stage_edge_lim < min(stage_centroid_values[k], zl +hc)){
-                //    stage_edge_lim = min(stage_centroid_values[k], zl+hc);
-                //}
-                //bedslope_work = g*hc*(stage_edge_lim)*length-0.5*g*max(stage_edge_lim-zl,0.)*(stage_edge_lim-zl)*length;
-
-                // Bedslope approx 3
-                //bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
-                //
-            //}else{
-            //   // Treat nearly dry cells 
-            //   bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
-            //   //
-            //   //bedslope_work = -pressure_flux*length;
-            //   xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
-            //   ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
-            //   //xmom_explicit_update[k] = 0.0;
-            //   //ymom_explicit_update[k] = 0.0;
-
-            ////}
-
             already_computed_flux[ki] = call; // #k Done
 
@@ -511 +406 @@
             // Update neighbour n with same flux but reversed sign
             if (n >= 0) {
-                stage_explicit_update[n] += edgeflux[0];
-                xmom_explicit_update[n] += edgeflux[1];
-                ymom_explicit_update[n] += edgeflux[2];
-                //Add bed slope term here
-                //if(hc_n>-9999.0){
-                //if(stage_centroid_values[n] > max_bed_edgevalue[n]){
-                    // Bedslope approx 1:
-                if(hc_n>H0){
-                    //bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
-                    bedslope_work = length*g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux*length;
-                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
-                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
-                }else{
-                    //bedslope_work = 0.;//-pressure_flux*length;//g*length*(0.*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
-                    xmom_explicit_update[n] *= 0.;
-                    ymom_explicit_update[n] *= 0.;
-                }
-                //if(k==229| n==229){
-                //    printf("-n-pres,xmom,bedwrk, %20.20e, %20.20e, %20.20e \n", pressure_flux*length, xmom_explicit_update[n], bedslope_work);
-                //}
-
-                    //
-                    // Bedslope approx 2:
-                    //stage_edge_lim = qr[0];
-                    //if(stage_edge_lim > max(stage_centroid_values[n], zr +hc_n)){
-                    //    stage_edge_lim = max(stage_centroid_values[n], zr+hc_n);
-                    //}
-                    //if(stage_edge_lim < min(stage_centroid_values[n], zr +hc_n)){
-                    //    stage_edge_lim = min(stage_centroid_values[n], zr+hc_n);
-                    //}
-                    //bedslope_work = g*hc_n*(stage_edge_lim)*length-0.5*g*max(stage_edge_lim-zr,0.)*(stage_edge_lim-zr)*length;
-                    //
-                    // Bedslope approx 3
-                    //bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
-                    //
-                //}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;
-                //    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
-                //    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
-
-                //    //xmom_explicit_update[n] = 0.0;
-                //    //ymom_explicit_update[n] = 0.0;
-                //}
+
+                edgeflux_store[nm3 + 0 ] = edgeflux[0];
+                edgeflux_store[nm3 + 1 ] = edgeflux[1];
+                edgeflux_store[nm3 + 2 ] = edgeflux[2];
+
+                bedslope_work = length*(g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux);
+                //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0];  
+                pressuregrad_store[nm] = bedslope_work;
+
 
                 already_computed_flux[nm] = call; // #n Done
@@ -561 +420 @@
 
             // Update timestep based on edge i and possibly neighbour n
-            if (tri_full_flag[k] == 1) {
+            // GD HACK:
+            // FIXME: Presently this 'hacked' to only update the timestep on
+            // every 2nd call.  This works for rk2 timestepping only, and is
+            // needed for correct wetting and drying treatment
+            if ((tri_full_flag[k] == 1) && (call%2==0)) {
                 if (max_speed > epsilon) {
                     // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
@@ -573 +436 @@
                     }
 
-                    // Ted Rigby's suggested less conservative version
-                    //if(n>=0){
-                    //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
-                    //}else{
-                    //  timestep = min(timestep, radii[k]/max_speed);
-                    //}
                 }
             }
@@ -584 +441 @@
         } // End edge i (and neighbour n)
 
-        //if(k==2082){
-        //    printf("xmom_explicit: %20.20e \n", xmom_explicit_update[k]);
-        //}
+        // Keep track of maximal speeds
+        max_speed_array[k] = max_speed;
+
+    } // End triangle k
+ 
+    // GD HACK 
+    // Limit edgefluxes, for mass conservation near wet/dry cells
+    for(k=0; k< number_of_elements; k++){
+        // Only check cells that are **near dry**
+        // Found I could get problems by checking all cells
+        if(stage_centroid_values[k] <= max_bed_edgevalue[k] + H0){
+            hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.);
+            // Loop over every edge
+            for(i = 0; i<3; i++){
+                if(i==0){
+                    // Add up the outgoing flux through the cell
+                    outgoing_flux_edges=0.0;
+                    for(useint=0; useint<3; useint++){
+                        if(edgeflux_store[3*(3*k+useint)]< 0.){
+                            //outgoing_flux_edges+=1.0;
+                            outgoing_flux_edges+=(edgeflux_store[3*(3*k+useint)]*timestep);
+                        }
+                    }
+                }
+
+                ki=3*k+i;   
+                ki2=ki*2;
+                ki3 = ki*3;
+                
+                // Prevent outflow from 'seriously' dry cells
+                // Idea: The cell will not go dry if:
+                // total_outgoing_flux <= cell volume = Area_triangle*hc
+                if((edgeflux_store[ki3]< 0.0) && (outgoing_flux_edges<0.)){
+                    
+                    // This bound could be improved (e.g. we could actually sum
+                    // the + and - fluxes and check if they are too large).
+                    // However, the advantage is that we don't have to worry
+                    // about subsequent changes to the + edgeflux caused by
+                    // constraints associated with neighbouring triangles.
+                    tmp = (areas[k]*hc)/(fabs(outgoing_flux_edges)+1.0e-10) ;
+                    if(tmp< 1.0){
+                        // Idea: we should be close to conserving mass if this holds, 
+                        // --- consider CFL condition. 
+                        tmp2 = min(hc/((stage_edge_values[ki]-bed_edge_values[ki])+1.0e-10), 1.0);
+                        //tmp2=(xmom_centroid_values[k]*xmom_centroid_values[k]+ymom_centroid_values[k]*ymom_centroid_values[k]);
+                        //tmp2/=(xmom_edge_values[ki]*xmom_edge_values[ki]+ ymom_edge_values[ki]*ymom_edge_values[ki]+1.0e-10);
+                        //tmp2 = sqrt(tmp2);
+                        //tmp2 = min(tmp2, 1.0);
+                        edgeflux_store[ki3+0]*=tmp2;
+                        edgeflux_store[ki3+1]*=tmp2;
+                        edgeflux_store[ki3+2]*=tmp2;
+
+                        // Compute neighbour edge index
+                        n = neighbours[ki];
+                        if(n>=0){
+                            nm = 3*n + neighbour_edges[ki];
+                            nm3 = nm*3;
+                            edgeflux_store[nm3+0]*=tmp2;
+                            edgeflux_store[nm3+1]*=tmp2;
+                            edgeflux_store[nm3+2]*=tmp2;
+                        }
+                    }
+                }
+            }
+        }
+     }
+
+
+    // Now add up stage, xmom, ymom explicit updates
+    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++){
+            ki=3*k+i;   
+            ki2=ki*2;
+            ki3 = ki*3;
+
+            // GD HACK
+            // Option to limit advective fluxes
+            if(hc > H0){
+                stage_explicit_update[k] += edgeflux_store[ki3+0];
+                // Cut these terms for shallow flows, and protect stationary states from round-off error
+                xmom_explicit_update[k] += edgeflux_store[ki3+1];
+                ymom_explicit_update[k] += edgeflux_store[ki3+2];
+            }else{
+                stage_explicit_update[k] += edgeflux_store[ki3+0];
+            }
+
+
+            // GD HACK
+            // Compute bed slope term
+            if(hc>H0){
+                xmom_explicit_update[k] -= normals[ki2]*pressuregrad_store[ki];
+                ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_store[ki];
+            }else{
+                xmom_explicit_update[k] *= 0.;
+                ymom_explicit_update[k] *= 0.;
+            }
+
+        } // end edge i
 
         // Normalise triangle k by area and store for when all conserved
@@ -594 +551 @@
         xmom_explicit_update[k] *= inv_area;
         ymom_explicit_update[k] *= inv_area;
-
-        //if(k==2082){
-        //    printf(" inv_area: %20.20e \n", inv_area);
-        //    printf("xmom_explicit: %20.20e \n", xmom_explicit_update[k]);
-        //}
-        // Keep track of maximal speeds
-        max_speed_array[k] = max_speed;
-
-    } // End triangle k
-  
-     
-    //bedtop=0.;
+    }  // end cell k
+
+    // Look for 'dry' edges, and set momentum (& component of momentum_update)
+    // in that direction to zero. If we don't do this, then it is possible for
+    // the bedslope term normal to that edge to be nonzero, and for momentum to
+    // keep 'building' up in the cell without being able to flow out.
+    //
     //for(k=0; k<number_of_elements; k++){
-    //      bedtop = bedtop + stage_explicit_update[k]*areas[k];
+    //    // Ignore "Dry" centroid cells, which will automatically have their
+    //    // momentum set to zero (in _protect)
+    //    if(stage_centroid_values[k] - bed_centroid_values[k] > H0){
+    //        for(i=0;i<3;i++){
+    //            ki=3*k+i;
+    //            if(stage_edge_values[ki] - bed_edge_values[ki] <=0.){
+    //              // Compute magnitude of momentum_update in direction normal to edge
+    //              tmp =  xmom_explicit_update[k]*normals[2*ki] + ymom_explicit_update[k]*normals[2*ki+1];
+    //              if(tmp > 0.){
+    //                  // Set component of momentum_update normal to edge to zero
+    //                  // [equivalently, subtract component of momentum_update normal to edge]
+    //                  xmom_explicit_update[k] -= tmp*normals[2*ki]; 
+    //                  ymom_explicit_update[k] -= tmp*normals[2*ki+1];  
+    //                  
+    //              }
+    //              
+    //              // Compute magnitude of momentum in direction normal to edge
+    //              //tmp =  xmom_centroid_values[k]*normals[2*ki] + ymom_centroid_values[k]*normals[2*ki+1];
+    //              //if(tmp > 0.){
+    //              //    // Set component of momentum normal to edge to zero
+    //              //    // [equivalently, subtract component of momentum normal to edge]
+    //              //    xmom_centroid_values[k] -= tmp*normals[2*ki]; 
+    //              //    ymom_centroid_values[k] -= tmp*normals[2*ki+1];  
+    //              //    
+    //              //}
+    //           } 
+    //        }
+    //    }
     //}
-    //printf("Bedtop = %20.20e \n ", bedtop);
-    //    for(i=0;i<number_of_elements;i++){
-    //        bedtop=max(bedtop, fabs(ymom_edge_values[3*k+i]));
-    //    }
+        
     //   
     //    
@@ -634 +610 @@
     free(max_bed_edgevalue);
     //free(min_bed_edgevalue);
+    free(edgeflux_store);
+    free(pressuregrad_store);
 
     return timestep;
@@ -668 +646 @@
       hc = wc[k] - zc[k];
       // Definine the maximum bed edge value on triangle k.
-      bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
-
-      if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
-        // Set momentum to zero and ensure h is non negative
+      //bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
+
+      //if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
+      //  xmomc[k]*=0.1;
+      //  ymomc[k]*=0.1;
+      //}
+
+
+      //if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
+      if (hc < minimum_allowed_height) {
+            // Set momentum to zero and ensure h is non negative
             xmomc[k] = 0.0;
             ymomc[k] = 0.0;
@@ -689 +674 @@
 
                  wc[k] = max(wc[k], bmin); 
-                 printf(" mass_add, %f, %f, %f,%d\n", xc[k], yc[k], wc[k], k) ;
+                 printf(" mass_add, %f, %f, %f,%f,%f,%d\n", xc[k], yc[k], wc[k],zc[k],wc[k]-zc[k], k) ;
              
 
@@ -804 +789 @@
   double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
   double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
-  double hc, h0, h1, h2, beta_tmp, hfactor, xtmp, ytmp, weight;
+  double hc, h0, h1, h2, beta_tmp, hfactor, xtmp, ytmp, weight, tmp;
   double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vect_norm, l1, l2;
   
@@ -844 +829 @@
       }
 
+  //
   // Compute some useful statistics on wetness/dryness
+  //
   for(k=0; k<number_of_elements;k++){ 
-      count_wet_neighbours[k]=0;
       cell_wetness_scale[k] = 0.;
       // Check if cell k is wet
       if(stage_centroid_values[k] > elevation_centroid_values[k] + minimum_allowed_height+epsilon){
+      //if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
           cell_wetness_scale[k] = 1.;  
       }
+  }
+  //
+  for(k=0; k<number_of_elements;k++){ 
       // Count the wet neighbours
+      count_wet_neighbours[k]=0;
       for (i=0; i<3; i++){
         ktmp = surrogate_neighbours[3*k+i];              
-        if(stage_centroid_values[ktmp] > max_elevation_edgevalue[ktmp]){
+        if( (ktmp >= 0) && cell_wetness_scale[ktmp]>0.){
+            count_wet_neighbours[k]+=1;
+        }else if(ktmp<=0){
+            // Boundary condition -- FIXME: assume wet
             count_wet_neighbours[k]+=1;
         }
+         
         
       }
@@ -865 +860 @@
   // First treat all 'fully wet' cells, then treat 'partially dry' cells
   for(k_wetdry=0; k_wetdry<2; k_wetdry++){
-      //if(((stage_centroid_values[k] > max_elevation_edgevalue[k]))==k_wetdry){
-      //if(((stage_centroid_values[k] < elevation_centroid_values[k]+minimum_allowed_height))==k_wetdry){
-      // For partially wet cells, we now know that the edges of neighbouring
-      // fully wet cells have been defined
 
       // Begin extrapolation routine
@@ -935 +926 @@
           //==============================================
           // Number of boundaries <= 1
+          // 'Typical case'
           //==============================================    
         
@@ -953 +945 @@
           //if(stage_centroid_values[k2]<=max_elevation_edgevalue[k2]+epsilon){
           //if(stage_centroid_values[k2]<=elevation_centroid_values[k2]+minimum_allowed_height+epsilon){
-          if(cell_wetness_scale[k2]==0.0){
+          if(k2<0 || cell_wetness_scale[k2]==0.0){
               k2 = k ;
           }
           //if(stage_centroid_values[k0]<=max_elevation_edgevalue[k0]+epsilon){
           //if(stage_centroid_values[k0]<=elevation_centroid_values[k0]+minimum_allowed_height+epsilon){
-          if(cell_wetness_scale[k0]==0.0){
+          if(k0 < 0 || cell_wetness_scale[k0]==0.0){
               k0 = k ;
           }
           //if(stage_centroid_values[k1]<=max_elevation_edgevalue[k1]+epsilon){
           //if(stage_centroid_values[k1]<=elevation_centroid_values[k1]+minimum_allowed_height+epsilon){
-          if(cell_wetness_scale[k1]==0.0){
+          if(k1 < 0 || cell_wetness_scale[k1]==0.0){
               k1 = k ;
           }
@@ -1003 +995 @@
           area2 = dy2*dx1 - dy1*dx2; 
            
-          // Treat triangles with zero or 1 wet neighbours. 
-          if ((area2 <= 0)|(cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
+          // Treat triangles with zero or 1 wet neighbours (area2 <=0.) 
+          if ((area2 <= 0.))// || (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
           {
             //if((cell_wetness_scale[k]==1.0) | (k==k0)&&(k==k1)&&(k==k2)){
-              stage_edge_values[k3]   = stage_centroid_values[k];
-              stage_edge_values[k3+1] = stage_centroid_values[k];
-              stage_edge_values[k3+2] = stage_centroid_values[k];
+            //if(count_wet_neighbours[k] > 5.){
+                stage_edge_values[k3]   = stage_centroid_values[k];
+                stage_edge_values[k3+1] = stage_centroid_values[k];
+                stage_edge_values[k3+2] = stage_centroid_values[k];
+            //}else{
             //}else{
               // Dry cell with a single 'wet' neighbour
@@ -1028 +1022 @@
               //if( (k==k0) & (k==k1) & (k==k2)){
               //    // Isolated wet cell  -- use constant depth extrapolation for stage
-              //    //stage_edge_values[k3]   = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3];
-              //    //stage_edge_values[k3+1] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1];
-              //    //stage_edge_values[k3+2] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2];
-
+              //    stage_edge_values[k3]   = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3];
+              //    stage_edge_values[k3+1] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1];
+              //    stage_edge_values[k3+2] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2];
+              //}
               //    // Isolated wet cell -- constant depth extrapolation, but don't let water flow uphill
               //    //stage_edge_values[k3]   = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3]);
@@ -1122 +1116 @@
           
           // Calculate heights of neighbouring cells
-          hc = stage_centroid_values[k]  - elevation_centroid_values[k];
-          h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
-          h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
-          h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
-          hmin = min(min(h0, min(h1, h2)), hc);
-
-          hfactor = 0.0;
-          //if (hmin > 0.001) 
-          if (hmin > 0.) 
-          //if (hc>0.0)
-          {
+          //hc = stage_centroid_values[k]  - elevation_centroid_values[k];
+          //h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
+          //h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
+          //h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
+          //hmin = min(min(h0, min(h1, h2)), hc);
+          //// GD FIXME: No longer needed? 
+          //hfactor = 0.0;
+          ////if (hmin > 0.001) 
+          //if (hmin > 0.) 
+          ////if (hc>0.0)
+          //{
             hfactor = 1.0 ;//hmin/(hmin + 0.004);
             //hfactor=hmin/(hmin + 0.004);
-          }
+          //}
           
           //-----------------------------------
@@ -1156 +1150 @@
           b = dx1*dq2 - dx2*dq1;
           b *= inv_area2;
-          
+          // d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth
+          //hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0);
+          //tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 
+          //           ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03;
+          //tmp /=max(hc*hc*hc, 1.0e-09);
+          //a = -xmom_centroid_values[k]*tmp;
+          //b = -ymom_centroid_values[k]*tmp;
           // Calculate provisional jumps in stage from the centroid 
           // of triangle k to its vertices, to be limited
@@ -1167 +1167 @@
           // from the centroid to the min and max
           find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
-          
+        
           beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
         
@@ -1173 +1173 @@
           // Limit the gradient
           // This is more complex for stage than other quantities
-          if((count_wet_neighbours[k]>0)){//&&(k_wetdry==0)){
+          //if((count_wet_neighbours[k]>0)&&(cell_wetness_scale[k]=1.0)){//&&(k_wetdry==0)){
+          //if(stage_centroid_values[k] > max_elevation_edgevalue[k]){
               limit_gradient(dqv, qmin, qmax, beta_tmp);
               stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
               stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
               stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
-          }else{
-          ////if(count_wet_neighbours[k]==0){
-          //    //weight=max(cell_wetness_scale[k] - (stage_centroid_values[k]-elevation_centroid_values[k]), 0.);
-          //    //weight/=(max_elevation_edgevalue[k] - elevation_centroid_values[k]);
-          //    //weight=min(weight, 1.0);
-          //    //weight==0;
-          //    //weight=0.;
-          //    // 'Shallow flow on steep slope'
-          //    //limit_gradient(dqv, qmin, qmax, beta_tmp);
-              stage_edge_values[k3+0] =stage_centroid_values[k] + dqv[0];
-              stage_edge_values[k3+1] =stage_centroid_values[k] + dqv[1];
-              stage_edge_values[k3+2] =stage_centroid_values[k] + dqv[2];
-
-          //    // There exists a neighbouring bed cell with elevation > minimum
-          //    // neighbouring stage value
-
-          //    // We have to be careful in this situation. Limiting the stage gradient can
-          //    // cause problems for shallow flows down steep slopes (rainfall-runoff type problems)
-          //    // Previously 'balance_deep_and_shallow' was used to deal with this issue
-          //    //for(i=0; i<3; i++){
-          //    //    stage_edge_values[k3+i] = weight*stage_edge_values[k3+i] + 
-          //    //                             (1.0-weight)*(stage_centroid_values[k] -elevation_centroid_values[k]
-          //    //                                           +elevation_edge_values[k3+i]);
-          //    //}
-          }
+          // IDEA: extrapolate assuming slope is that of steady flow??
+          // GD HACK - FIXME
+          // Note that 'near-dry' cells which don't exchange mass with
+          // neighbouring cells have to revert to having a constant stage
+          // extrapolation, or else they will have a residual bed slope term.
+          // This helps! -- But probably makes rainfall-runoff problems worse.
+          //ii = (stage_edge_values[k3+0] > elevation_edge_values[k3+0])+ 
+          //     (stage_edge_values[k3+1] > elevation_edge_values[k3+1])+ 
+          //     (stage_edge_values[k3+2] > elevation_edge_values[k3+2]);
+          //if(ii<2){
+          //    stage_edge_values[k3+0] = stage_centroid_values[k] ;
+          //    stage_edge_values[k3+1] = stage_centroid_values[k] ;
+          //    stage_edge_values[k3+2] = stage_centroid_values[k] ;
+          //}
+          //}else{
+          // d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth
+          //    hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0);
+          //    tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 
+          //               ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03;
+          //    tmp /=max(hc*hc*hc, 1.0e-09);
+          //    a = -xmom_centroid_values[k]*tmp;
+          //    b = -ymom_centroid_values[k]*tmp;
+          //    dqv[0] = a*dxv0 + b*dyv0;
+          //    dqv[1] = a*dxv1 + b*dyv1;
+          //    dqv[2] = a*dxv2 + b*dyv2;
+          //////if(count_wet_neighbours[k]==0){
+          ////    //weight=max(cell_wetness_scale[k] - (stage_centroid_values[k]-elevation_centroid_values[k]), 0.);
+          ////    //weight/=(max_elevation_edgevalue[k] - elevation_centroid_values[k]);
+          ////    //weight=min(weight, 1.0);
+          ////    //weight==0;
+          //      weight=0.;
+          ////    // 'Shallow flow on steep slope'
+              //limit_gradient(dqv, qmin, qmax, beta_tmp);
+              //stage_edge_values[k3+0] =stage_centroid_values[k] + dqv[0];
+              //stage_edge_values[k3+1] =stage_centroid_values[k] + dqv[1];
+              //stage_edge_values[k3+2] =stage_centroid_values[k] + dqv[2];
+          //
+          ////    // There exists a neighbouring bed cell with elevation > minimum
+          ////    // neighbouring stage value
+
+          ////    // We have to be careful in this situation. Limiting the stage gradient can
+          ////    // cause problems for shallow flows down steep slopes (rainfall-runoff type problems)
+          ////    // Previously 'balance_deep_and_shallow' was used to deal with this issue
+          //    for(i=0; i<3; i++){
+          //        stage_edge_values[k3+i] = weight*stage_edge_values[k3+i] + 
+          //                                 (1.0-weight)*(stage_centroid_values[k] -elevation_centroid_values[k]
+          //                                               +elevation_edge_values[k3+i]);
+          //    }
+          //}
            
           //-----------------------------------
@@ -1240 +1265 @@
           //if((k!=k0)&(k!=k1)&(k!=k2))
           //if(stagemin>=bedmax){
-          if((count_wet_neighbours[k]>0) && 
-             (cell_wetness_scale[k]==1.0)){
-            //(stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
+          if((count_wet_neighbours[k]>0)){ // && 
+            // (cell_wetness_scale[k]==1.0)){
+          //if(stage_centroid_values[k]>(minimum_allowed_height+max_elevation_edgevalue[k])){
             limit_gradient(dqv, qmin, qmax, beta_tmp);
           }else{
@@ -1293 +1318 @@
           // Limit the gradient
           //if(stagemin>=bedmax){
-          if((count_wet_neighbours[k]>0)&&
-             (cell_wetness_scale[k]==1.0)){
+          if((count_wet_neighbours[k]>0)){//&&
+             //(cell_wetness_scale[k]==1.0)){
             //(stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
+          //if(stage_centroid_values[k]>(minimum_allowed_height+max_elevation_edgevalue[k])){
             limit_gradient(dqv, qmin, qmax, beta_tmp);
           }else{
@@ -1635 +1661 @@
     *max_speed_array, //Keeps track of max speeds for each triangle
     *stage_centroid_values,
+    *xmom_centroid_values,
+    *ymom_centroid_values,
     *bed_centroid_values,
     *bed_vertex_values;
@@ -1642 +1670 @@
     
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOiOOO",
+  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOiOOOOO",
             &timestep,
             &epsilon,
@@ -1667 +1695 @@
             &optimise_dry_cells,
             &stage_centroid_values,
+            &xmom_centroid_values,
+            &ymom_centroid_values,
             &bed_centroid_values,
             &bed_vertex_values)) {
@@ -1695 +1725 @@
   CHECK_C_CONTIG(max_speed_array);
   CHECK_C_CONTIG(stage_centroid_values);
+  CHECK_C_CONTIG(xmom_centroid_values);
+  CHECK_C_CONTIG(ymom_centroid_values);
   CHECK_C_CONTIG(bed_centroid_values);
   CHECK_C_CONTIG(bed_vertex_values);
@@ -1729 +1761 @@
                      optimise_dry_cells, 
                      (double*) stage_centroid_values -> data, 
+                     (double*) xmom_centroid_values -> data, 
+                     (double*) ymom_centroid_values -> data, 
                      (double*) bed_centroid_values -> data,
                      (double*) bed_vertex_values -> data);

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

@@ -14 +14 @@
 #from anuga import *
 #from balanced_basic import *
-#from balanced_dev import *
-from anuga_tsunami import *
+from balanced_dev import *
+#from anuga_tsunami import *
 #from balanced_basic.swb2_domain import *
 #from balanced_basic.swb2_domain import Domain

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

@@ -51 +51 @@
 # Setup Algorithm
 #------------------------------------------------------------------------------
-domain.set_timestepping_method('rk2')
-domain.set_default_order(2)
+#domain.set_timestepping_method('rk2')
+#domain.set_default_order(2)
 
 print domain.get_timestepping_method()
 
-domain.use_edge_limiter = True
+#domain.use_edge_limiter = True
 #domain.use_edge_limiter = False
-domain.tight_slope_limiters = True
-domain.use_centroid_velocities = False
+#domain.tight_slope_limiters = True
+#domain.use_centroid_velocities = False
 
-domain.CFL = 1.0
+#domain.CFL = 1.0
 
 #domain.beta_w      = 0.6

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_20120831_172309/dam_break.sww', 0.001)
+p_dev = util.get_output('dam_break_20130219_202649/dam_break.sww', 0.001)
 p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
 

trunk/anuga_work/development/gareth/tests/okushiri/run_okushiri.py

@@ -24 +24 @@
 import anuga
 import project
-#from balanced_dev import *
-from anuga_tsunami import *
+from balanced_dev import *
+#from anuga_tsunami import *
 
 def main(elevation_in_mesh=False):

trunk/anuga_work/development/gareth/tests/parabolic/parabolaplot.py

@@ -11 +11 @@
 import util
 
-p=util.get_output('parabola_v2.sww', 0.01)
+p=util.get_output('parabola_v2.sww', 0.001)
 p2=util.get_centroids(p, velocity_extrapolation=True)
 

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

@@ -56 +56 @@
 Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
 Bd=anuga.Dirichlet_boundary([-0.2,0.,0.])       # Constant boundary values -- not used in this example
-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.
+#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.
 
 #----------------------------------------------

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

@@ -14 +14 @@
 #from swb2_domain import *
 #from balanced_basic import *
-#from balanced_dev import *
-from anuga_tsunami import *
+from balanced_dev import *
+#from anuga_tsunami import *
+
 #---------
 #Setup computational domain
@@ -23 +24 @@
 domain=Domain(points,vertices,boundary)    # Create Domain
 domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
-domain.set_datadir('.')                          # Use current folder
-domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
-#domain.set_store_vertices_uniquely(True)
+
 #------------------
 # Define topography
 #------------------
 scale_me=1.0
-
 #domain.minimum_allowed_height=domain.minimum_allowed_height*scale_me # Seems needed to make the algorithms behave
 
@@ -37 +35 @@
 
 def stagefun(x,y):
-    stge=-0.2*scale_me #+0.01*(x>0.9) 
+    stge=-0.2*scale_me #+0.1*(x>0.9) 
     #topo=topography(x,y) 
     return stge#*(stge>topo) + (topo)*(stge<=topo)
@@ -53 +51 @@
 domain.get_quantity('stage').smooth_vertex_values()
 
+#print domain.forcing_terms
+#assert 0==1
 # Experiment with rain.
 # rainin = anuga.shallow_water.forcing.Rainfall(domain, rate=0.001) #, center=(0.,0.), radius=1000. )
@@ -64 +64 @@
 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
+    return -0.1 #(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,
@@ -72 +72 @@
 # Associate boundary tags with boundary objects
 #----------------------------------------------
-domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
+domain.set_boundary({'left': Br, 'right': Bt2, 'top': Br, 'bottom':Br})
 
 #------------------------------
@@ -78 +78 @@
 #------------------------------
 
-for t in domain.evolve(yieldstep=0.2,finaltime=40.0):
+for t in domain.evolve(yieldstep=0.200,finaltime=40.00):
     print domain.timestepping_statistics()
     xx = domain.quantities['xmomentum'].centroid_values
@@ -84 +84 @@
     dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values 
     dd_raw=1.0*dd
-    dd = (dd)*(dd>1.0e-03)+1.0e-03
+    dd = dd*(dd>1.0e-03)+1.0e-03*(dd<=1.0e-03)
     vv = ( (xx/dd)**2 + (yy/dd)**2)**0.5
     vv = vv*(dd>1.0e-03)

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

@@ -15 +15 @@
 #------------------
 #v=(p.y==0.5)
-v=(p2.y==p2.y[80])
+myind=(abs(p2.y-0.5)).argmin()
+v=(p2.y==p2.y[myind])
 
 #--------------------

trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py

@@ -29 +29 @@
 # Setup computational domain
 #------------------------------------------------------------------------------
-points, vertices, boundary = rectangular_cross(10, 10,
-                                len1=100.0, len2=100.0) # Mesh
+points, vertices, boundary = rectangular_cross(14, 10,
+                                len1=140.0, len2=100.0) # Mesh
 domain = Domain(points, vertices, boundary) # Create domain
 domain.set_name('channel_SU_2_v2') # Output name
-#domain.set_store_vertices_uniquely(True)
-#domain.CFL=1.0
 #------------------------------------------------------------------------------
 # Setup initial conditions
 #------------------------------------------------------------------------------
 def topography(x, y):
-        return -x/10. #+ numpy.sin(x/10.) +abs(y-50.)/100. # linear bed slope
+        return -x/10. + 1.*(numpy.sin(x/10.) +abs(y-50.)/10.) -0.*(x>80.) # linear bed slope
 
 def stagetopo(x,y):
-    stg= -x/10. +0.06309625 # Constant depth
+    stg= topography(x,y) + 0.163 #-x/10. +0.06309625 -50.*(x>80.)# Constant depth
     #topo=topography(x,y)
     return stg#*(stg>topo) + topo*(stg<=topo)
@@ -68 +66 @@
 #Bo = anuga.Transmissive_boundary(domain)
 #Bd = anuga.Dirichlet_boundary([-20., 0., 0.])
-domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})