"""Test script that shows issue with outflow forcing in very shallow water.
There is currently no way of knowing how much water was actually removed

"""
                    
from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
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.
width = 10.

dx = dy = 2  # 1 or 2 OK for this test
points, vertices, boundary = rectangular_cross(int(length/dx),
                                               int(width/dy),
                                               len1=length, 
                                               len2=width)
domain = Domain(points, vertices, boundary)   
domain.set_name('test_outflow_conservation')  # Output name
domain.set_default_order(2)
        

# Flat surface
stage = 0.25  # Values greater than about 0.5 work. Values at 0.3 or less fail 

domain.set_quantity('elevation', 0)
domain.set_quantity('stage', stage)
domain.set_quantity('friction', 0)

Br = Reflective_boundary(domain)
domain.set_boundary({'left': Br, 'right': Br, 'bottom': Br, 'top': Br})

# 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(I)
initial_volume = domain.quantities['stage'].get_integral()
predicted_volume = initial_volume
dt = 0.05
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 + I.rate/pi/100/dt