Changeset 7345 for anuga_core/documentation/user_manual/demos/channel_variable.py

Timestamp:

Aug 9, 2009, 5:13:36 PM (15 years ago)

Author:

ole

Message:

Working variable bed example - storing time dependent elevation has not been attempted yet.

File:

• anuga_core/documentation/user_manual/demos/channel_variable.py (modified) (5 diffs)

anuga_core/documentation/user_manual/demos/channel_variable.py

@@ -32 +32 @@
     """Complex topography defined by a function of vectors x and y."""
 
-    z = -x/10
+    z = -x/100
 
     N = len(x)
@@ -41 +41 @@
 
         # Constriction
-        if 27 < x[i] < 29 and y[i] > 3:
-            z[i] += 2
+        #if 27 < x[i] < 29 and y[i] > 3:
+        #    z[i] += 2
 
         # Pole
@@ -50 +50 @@
     return z
 
-def pole_increments(x,y):
+def pole_increment(x,y):
     """This provides a small increment to a pole located mid stream
     For use with variable elevation data
@@ -60 +60 @@
     for i in range(N):
         # Pole
-        if (x[i] - 34)**2 + (y[i] - 2)**2 < 0.4**2:
-            z[i] += 0.01
-
+        if (x[i] - 20)**2 + (y[i] - 2)**2 < 0.4**2:
+            z[i] += 0.05
     return z
     
@@ -82 +81 @@
 # Evolve system through time
 #------------------------------------------------------------------------------
-for t in domain.evolve(yieldstep=0.1, finaltime=16.0):
+
+growing = False
+shrinking = False
+done = False
+for t in domain.evolve(yieldstep=0.1, finaltime=30.0):
     print domain.timestepping_statistics()
 
-    if domain.get_quantity('stage').\
-           get_values(interpolation_points=[[10, 2.5]]) > 0:
-        print 'Stage > 0: Changing elevation data'
+    #w = domain.get_quantity('stage').\
+    #    get_values(interpolation_points=[[18, 2.5]])
+    #print 'Level at gauge point = %.2fm' % w
+           
+    z = domain.get_quantity('elevation').\
+        get_values(interpolation_points=[[20, 2]])           
+    print 'Elevation at pole location = %.2fm' % z           
+
+    # Start variable elevation after 10 seconds    
+    if t > 10 and not (shrinking or growing or done):
+        growing = True
+           
+    # Turn around when pole has reached a height of 2 m
+    if z >= 2 and growing:
+        growing = False
+        shrinking = True
         
+    # Stop changing when pole has shrunk to 1 m
+    if z <= 1 and shrinking:
+        done = True
+        shrinking = growing = False
+
+    # Grow or shrink               
+    if growing:       
         domain.add_quantity('elevation', pole_increment)
+        
+    if shrinking:    
+        domain.add_quantity('elevation', lambda x,y: -pole_increment(x,y))    
+