Changeset 8355 for trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

Timestamp:

Mar 8, 2012, 5:07:31 PM (13 years ago)

Author:

davies

Message:

balanced_dev major bugfix

File:

• trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c (modified) (5 diffs)

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

@@ -470 +470 @@
             }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;
@@ -512 +512 @@
                 }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;
@@ -782 +782 @@
   double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
   double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
-  double dqv[3], qmin, qmax, hmin, hmax;
+  double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
   double hc, h0, h1, h2, beta_tmp, hfactor;
   double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2;
@@ -867 +867 @@
       k1 = surrogate_neighbours[k3 + 1];
       k2 = surrogate_neighbours[k3 + 2];
-      // If the mean stage value in neighbouring triangles is less than
-      // the bed edge value at the common edge, then don't use this triangle
-      // in the extrapolation. 
+     
+      // Test to see whether we accept the surrogate neighbours 
       // Note that if ki is replaced with k in more than 1 neighbour, then the
       // triangle area will be zero, and a first order extrapolation will be
       // used
-      if(0.5*(stage_centroid_values[k2]+stage_centroid_values[k])<=elevation_edge_values[k2]+minimum_allowed_height){
-          k2 = k ;
-      }
-      if(0.5*(stage_centroid_values[k0]+stage_centroid_values[k])<=elevation_edge_values[k0]+minimum_allowed_height){
-          k0 = k ;
-      }
-      if(0.5*(stage_centroid_values[k1]+stage_centroid_values[k])<=elevation_edge_values[k1]+minimum_allowed_height){
-          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;
+      } 
+     
       // Get the auxiliary triangle's vertex coordinates 
       // (really the centroids of neighbouring triangles)
@@ -916 +931 @@
           //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]   = 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];
+          stage_edge_values[k3+2] = stage_centroid_values[k];
           xmom_edge_values[k3]    = xmom_centroid_values[k];
           xmom_edge_values[k3+1]  = xmom_centroid_values[k];