Changeset 6142


Ignore:
Timestamp:
Jan 13, 2009, 10:51:33 AM (10 years ago)
Author:
ole
Message:

Developed protection against situation where requested culvert flows exceed
what is possible at the inlet. This works by computing the largest possible flow for each triangle and use that to reduce the overall flow if necessary.

Location:
anuga_core/source/anuga/culvert_flows
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/culvert_flows/culvert_class.py

    r6131 r6142  
    77
    88from anuga.utilities.numerical_tools import mean
    9 from anuga.utilities.numerical_tools import ensure_numeric
     9from anuga.utilities.numerical_tools import ensure_numeric, sign
    1010       
    1111from anuga.config import g, epsilon
     
    401401           
    402402
    403            
     403    def adjust_flow_for_available_water_at_inlet(self, Q, delta_t):
     404        """Adjust Q downwards depending on available water at inlet
     405        """
     406       
     407        # Short hands
     408        domain = self.domain       
     409        dq = domain.quantities               
     410        time = domain.get_time()
     411        I = self.inlet
     412
     413       
     414        # Find triangle with the smallest available flow
     415        max_Q = sys.maxint
     416        if delta_t > 0:
     417            for i in I.exchange_indices:
     418                stage = dq['stage'].get_values(location='centroids',
     419                                               indices=[i])[0]
     420                elevation = dq['elevation'].get_values(location='centroids',
     421                                                       indices=[i])[0]       
     422                depth = stage-elevation
     423                area = domain.areas[i]
     424
     425                # Possible rate based on this triangle
     426                Q_possible = depth*I.exchange_area/delta_t
     427                               
     428                # Use triangle with least depth to determine total flow
     429                max_Q = min(max_Q, Q_possible)
     430
     431           
     432        Q_reduced = sign(Q)*min(abs(Q), abs(max_Q))
     433       
     434        if abs(Q_reduced) < abs(Q):
     435            msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
     436            msg += 'greater than current volume (V) at inlet.\n'
     437            msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
     438            if self.verbose is True:
     439                print msg
     440            if hasattr(self, 'log_filename'):                   
     441                log_to_file(self.log_filename, msg)
     442       
     443        return Q_reduced   
     444                       
    404445           
    405446    def compute_rates(self, delta_t):
     
    544585       
    545586
    546         # Adjust Q downwards depending on available water at inlet
    547         # Experimental
    548         stage = self.inlet.get_quantity_values(quantity_name='stage')
    549         elevation = self.inlet.get_quantity_values(quantity_name='elevation')
    550         depth = stage-elevation
    551        
    552        
    553         V = 0
    554         for i, d in enumerate(depth):
    555             V += d * domain.areas[i]
    556        
    557         #Vsimple = mean(depth)*self.inlet.exchange_area # Current volume in exchange area 
    558         #print 'Q', Q, 'dt', delta_t, 'Q*dt', Q*delta_t, 'V', V, 'Vsimple', Vsimple
    559 
    560         dt = delta_t           
    561         if Q*dt > V:
    562        
    563             Q_reduced = 0.9*V/dt # Reduce with safety factor
    564            
    565             msg = '%.2fs: Computed extraction for this time interval (Q*dt) is ' % time
    566             msg += 'greater than current volume (V) at inlet.\n'
    567             msg += ' Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced)
    568            
    569             #print msg
    570            
    571             if self.verbose is True:
    572                 print msg
    573             if hasattr(self, 'log_filename'):                   
    574                 log_to_file(self.log_filename, msg)                                       
    575            
    576             Q = Q_reduced
     587        Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)
    577588       
    578589        self.inlet.rate = -Q
  • anuga_core/source/anuga/culvert_flows/test_culvert_class.py

    r6128 r6142  
    319319                               end_point0=[9.0, 2.5],
    320320                               end_point1=[13.0, 2.5],
    321                                trigger_depth=0.6, #FIXME: This causes test to pass, but smaller values do not
     321                               trigger_depth=0.01,
    322322                               verbose=False)
    323323
     
    346346            msg = 'Total volume has changed: Is %.2f m^3 should have been %.2f m^3'\
    347347                % (new_volume, ref_volume)
     348            if not allclose(new_volume, ref_volume):
     349                print msg
    348350            assert allclose(new_volume, ref_volume), msg
    349351   
Note: See TracChangeset for help on using the changeset viewer.