Changeset 8253

Timestamp:

Nov 21, 2011, 1:36:41 PM (13 years ago)

Author:

paul

Message:

sqpipe update_state now uses neighbouring states

Location:

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe

Files:

• sqpipe_domain.py (modified) (6 diffs)
• sqpipe_test_domain.py (modified) (1 diff)

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py

@@ -50 +50 @@
 class Sqpipe_domain(Generic_domain):
 
-    def __init__(self, coordinates, conserved_quantities, forcing_terms = [], boundary = None, state = None, update_state = True, bulk_modulus = 1.0, tagged_elements = None):
+    def __init__(self, coordinates, conserved_quantities, forcing_terms = [], boundary = None, state = None, update_state_flag = True, bulk_modulus = 1.0, tagged_elements = None):
 
         evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','top','stage']
@@ -70 +70 @@
             self.state = state            
 
-        self.update_state = update_state
+        self.update_state_flag = update_state_flag
 
         # Get computational parameters
@@ -126 +126 @@
             q = [q] # Turn argument into a list
 
-        #Check correcness    
+        #Check correctness    
         for quantity_name in q:
             msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
@@ -145 +145 @@
 
         # Update state before computing flux (if necessary)
-        if self.update_state :
-            h = self.quantities['height'].centroid_values
-            t = self.quantities['top'].centroid_values
-            self.state = numpy.where(h>t, 1, 0)
+        if self.update_state_flag:
+            self.update_state()
 
         # Import flux method and call it
@@ -154 +152 @@
         self.flux_timestep = compute_fluxes_ext(timestep,self)        
 
+    def update_state(self):
+        h = self.quantities['height'].centroid_values
+        t = self.quantities['top'].centroid_values
+
+        new_state = numpy.zeros_like(self.state)
+
+        # Update boundary state
+        new_state[0] = update_cell_state(0, [0, 1], h, t, self.state)
+        N = self.number_of_elements-1
+        new_state[N] = update_cell_state(N, [N-1, N], h, t, self.state)
+
+        # Update interior states
+        for i in range(1, N-1):
+            new_state[i] = update_cell_state(i, [i-1, i+1], h, t, self.state)
+
+        self.state = new_state
+        #self.state = numpy.where(h>=t, 1, 0)
 
     def distribute_to_vertices_and_edges(self):
@@ -171 +186 @@
             yield(t)
 
+# Auxillary methods
+
+# Method to update the state of a cell
+# cell is the index of the cell to update
+# neighbours is a list of the indices of it's neighbours
+# (for boundary cells, this could include the cell itself) 
+def update_cell_state(cell, neighbours, h, t, s):
+    # If height is bigger than top, then pressurise
+    if h[cell] >= t[cell]:
+        return 1
+    else:
+        # If the cell was pressurised, and all it's neighbours
+        # are pressurised, it remains pressurised even though
+        # height is less than top
+        # Otherwise, it's free surface
+        if all([s[i] for i in neighbours]):
+            return s[cell]
+        else:
+            return 0
+
 #=============== End of Shallow Water Domain ===============================

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_test_domain.py

@@ -17 +17 @@
                                state                = state,
                                bulk_modulus         = bulk_modulus,
-        #                       update_state         = False,
+        #                       update_state_flag         = False,
                                tagged_elements      = tagged_elements)