Changeset 6141


Ignore:
Timestamp:
Jan 13, 2009, 10:14:17 AM (10 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:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/debug/outflow_conservation_example.py

    r6122 r6141  
    77from anuga.shallow_water.shallow_water_domain import Domain, Inflow, Reflective_boundary
    88from math import pi, cos, sin
     9from Numeric import allclose
     10import sys
    911       
    1012length = 20.
     
    3335# Apply outflow
    3436rate = -2.0
     37I = Inflow(domain, rate=rate, center=(15,5), radius=1)       
     38I.requested_rate = I.rate
     39
    3540domain.forcing_terms = []       
    36 domain.forcing_terms.append(Inflow(domain, rate=rate, center=(15,5), radius=1))       
     41domain.forcing_terms.append(I)
    3742initial_volume = domain.quantities['stage'].get_integral()
    3843predicted_volume = initial_volume
     
    4045for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
    4146    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 = stage-elevation
     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   
    4281   
    4382    print t, 'Volume difference:', abs(volume-predicted_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
    4785           
Note: See TracChangeset for help on using the changeset viewer.