source: anuga_work/debug/outflow_conservation_example.py @ 6501

Last change on this file since 6501 was 6141, checked in by ole, 15 years ago

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

File size: 3.1 KB
Line 
1"""Test script that shows issue with outflow forcing in very shallow water.
2There is currently no way of knowing how much water was actually removed
3
4"""
5                   
6from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
7from anuga.shallow_water.shallow_water_domain import Domain, Inflow, Reflective_boundary
8from math import pi, cos, sin
9from Numeric import allclose
10import sys
11       
12length = 20.
13width = 10.
14
15dx = dy = 2  # 1 or 2 OK for this test
16points, vertices, boundary = rectangular_cross(int(length/dx),
17                                               int(width/dy),
18                                               len1=length, 
19                                               len2=width)
20domain = Domain(points, vertices, boundary)   
21domain.set_name('test_outflow_conservation')  # Output name
22domain.set_default_order(2)
23       
24
25# Flat surface
26stage = 0.25  # Values greater than about 0.5 work. Values at 0.3 or less fail
27
28domain.set_quantity('elevation', 0)
29domain.set_quantity('stage', stage)
30domain.set_quantity('friction', 0)
31
32Br = Reflective_boundary(domain)
33domain.set_boundary({'left': Br, 'right': Br, 'bottom': Br, 'top': Br})
34
35# Apply outflow
36rate = -2.0
37I = Inflow(domain, rate=rate, center=(15,5), radius=1)       
38I.requested_rate = I.rate
39
40domain.forcing_terms = []       
41domain.forcing_terms.append(I)
42initial_volume = domain.quantities['stage'].get_integral()
43predicted_volume = initial_volume
44dt = 0.05
45for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
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 = 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   
81   
82    print t, 'Volume difference:', abs(volume-predicted_volume) 
83    assert allclose (volume, predicted_volume)           
84    predicted_volume = predicted_volume + I.rate/pi/100/dt
85           
Note: See TracBrowser for help on using the repository browser.