Changeset 6141
 Timestamp:
 Jan 13, 2009, 10:14:17 AM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/debug/outflow_conservation_example.py
r6122 r6141 7 7 from anuga.shallow_water.shallow_water_domain import Domain, Inflow, Reflective_boundary 8 8 from math import pi, cos, sin 9 from Numeric import allclose 10 import sys 9 11 10 12 length = 20. … … 33 35 # Apply outflow 34 36 rate = 2.0 37 I = Inflow(domain, rate=rate, center=(15,5), radius=1) 38 I.requested_rate = I.rate 39 35 40 domain.forcing_terms = [] 36 domain.forcing_terms.append(I nflow(domain, rate=rate, center=(15,5), radius=1))41 domain.forcing_terms.append(I) 37 42 initial_volume = domain.quantities['stage'].get_integral() 38 43 predicted_volume = initial_volume … … 40 45 for t in domain.evolve(yieldstep = dt, finaltime = 5.0): 41 46 volume = domain.quantities['stage'].get_integral() 47 48 I.rate = I.requested_rate 49 50 delta_t = domain.timestep 51 if delta_t > 0: 52 max_Q = sys.maxint 53 for i in I.exchange_indices: 54 stage = domain.get_quantity('stage').get_values(location='centroids', 55 indices=[i])[0] 56 elevation = domain.get_quantity('elevation').get_values(location='centroids', 57 indices=[i])[0] 58 depth = stageelevation 59 area = domain.areas[i] 60 print 'exchange', I.exchange_area, 'timestep', delta_t 61 62 print i, 63 print ', area', area, ', depth', depth 64 print 'Requested rate for this triangle = %f m^3/s' % (I.rate/I.exchange_area*area) 65 66 Qi = depth*area/delta_t 67 print 'Possible rate for this triangle = %f m^3/s' % (Qi) 68 69 print 'Requested rate for exchange area = %f m^3/s' % (I.rate) 70 Qtot = Qi/area*I.exchange_area 71 #print 'Possible rate for exchange area = %f m^3/s' % (Qtot) 72 Qtot = 0.7*Qtot 73 74 max_Q = min(max_Q, Qtot) 75 76 77 78 I.rate = min(I.requested_rate, max_Q) 79 print 'Adjusted_rate for exchange_area = %f m^3/s' % I.rate 80 42 81 43 82 print t, 'Volume difference:', abs(volumepredicted_volume) 44 45 #assert allclose (volume, predicted_volume) 46 predicted_volume = predicted_volume + rate/pi/100/dt 83 assert allclose (volume, predicted_volume) 84 predicted_volume = predicted_volume + I.rate/pi/100/dt 47 85
Note: See TracChangeset
for help on using the changeset viewer.