Changeset 6141

Timestamp:

Jan 13, 2009, 10:14:17 AM (16 years ago)

Author:

ole

Message:

Tested out idea where Inflow rate is reduced by what is possibel on a triangle by triangle basis. This looks good!!

File:

• anuga_work/debug/outflow_conservation_example.py (modified) (3 diffs)

anuga_work/debug/outflow_conservation_example.py

@@ -7 +7 @@
 from anuga.shallow_water.shallow_water_domain import Domain, Inflow, Reflective_boundary
 from math import pi, cos, sin
+from Numeric import allclose
+import sys
         
 length = 20.
@@ -33 +35 @@
 # Apply outflow 
 rate = -2.0
+I = Inflow(domain, rate=rate, center=(15,5), radius=1)        
+I.requested_rate = I.rate
+
 domain.forcing_terms = []        
-domain.forcing_terms.append(Inflow(domain, rate=rate, center=(15,5), radius=1))        
+domain.forcing_terms.append(I)
 initial_volume = domain.quantities['stage'].get_integral()
 predicted_volume = initial_volume
@@ -40 +45 @@
 for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
     volume = domain.quantities['stage'].get_integral()
+
+    I.rate = I.requested_rate
+    
+    delta_t = domain.timestep
+    if delta_t > 0:
+        max_Q = sys.maxint
+        for i in I.exchange_indices:
+            stage = domain.get_quantity('stage').get_values(location='centroids', 
+                                                            indices=[i])[0]
+            elevation = domain.get_quantity('elevation').get_values(location='centroids', 
+                                                                    indices=[i])[0]        
+            depth = stage-elevation
+            area = domain.areas[i]
+            print 'exchange', I.exchange_area, 'timestep', delta_t
+            
+            print i,
+            print ', area', area, ', depth', depth 
+            print 'Requested rate for this triangle = %f m^3/s' % (I.rate/I.exchange_area*area)
+            
+            Qi = depth*area/delta_t
+            print 'Possible rate for this triangle = %f m^3/s' % (-Qi)
+            
+            print 'Requested rate for exchange area = %f m^3/s' % (I.rate)            
+            Qtot = Qi/area*I.exchange_area
+            #print 'Possible rate for exchange area = %f m^3/s' % (-Qtot)            
+            Qtot = 0.7*Qtot
+            
+            max_Q = min(max_Q, Qtot)
+
+            
+            
+        I.rate = -min(-I.requested_rate, max_Q)
+        print 'Adjusted_rate for exchange_area = %f m^3/s' % I.rate                    
+    
     
     print t, 'Volume difference:', abs(volume-predicted_volume) 
-    
-    #assert allclose (volume, predicted_volume)            
-    predicted_volume = predicted_volume + rate/pi/100/dt
+    assert allclose (volume, predicted_volume)            
+    predicted_volume = predicted_volume + I.rate/pi/100/dt