Changeset 7573 for anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

Timestamp:

Nov 27, 2009, 9:31:34 AM (13 years ago)

Author:

steve

Message:

Committing a version of shallow_water_balanced which passes it unit test
using a version of edge limiting which doesn't limit boundary edges. THis
is useful to allow linear functions to be reconstructed.

Had to play with the transmissive BC and use centroid values which is
set via domain set_centroid_transmissive_bc

File:

• anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -339 +339 @@
   int i, k, k2, k3, k6;
   long n;
-  double qmin, qmax, qn, qc;
+  double qmin, qmax, qn, qc, sign;
   double dq, dqa[3], phi, r;
   
@@ -361 +361 @@
     }
     
+    sign = 0.0;
+    if (qmin > 0) {
+      sign = 1.0;
+    } else if (qmax < 0) {
+      sign = -1.0;
+    }
+
     phi = 1.0;
-    for (i=0; i<3; i++) {    
-      r = 1.0;
-      
-      dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
-      dqa[i] = dq;                      //Save dq for use in updating vertex values
-      
-      if (dq > 0.0) r = (qmax - qc)/dq;
-      if (dq < 0.0) r = (qmin - qc)/dq;      
-      
-      
-      phi = min( min(r*beta, 1.0), phi);    
+    for (i=0; i<3; i++) {  
+      dq = edge_values[k3+i] - qc;      //Delta between edge and centroid values
+      dqa[i] = dq;                      //Save dq for use in updating vertex values  
+      
+      // FIXME SR 20091125 This caused problems in shallow_water_balanced
+      // commenting out problem
+      // Just limit non boundary edges so that we can reconstruct a linear function
+      if (neighbours[k3+i] >= 0) {
+        r = 1.0;
+      
+        if (dq > 0.0) r = (qmax - qc)/dq;
+        if (dq < 0.0) r = (qmin - qc)/dq;      
+            
+        phi = min( min(r*beta, 1.0), phi);
+        }
+
+      if (neighbours[k3+i] < 0) {
+        r = 1.0;
+      
+        if (dq > 0.0 && sign == -1.0 ) r = (0.0 - qc)/dq;
+        if (dq < 0.0 && sign ==  1.0 ) r = (0.0 - qc)/dq;      
+            
+        phi = min( min(r*beta, 1.0), phi);
+        }
+    
     }