1 | #!/usr/bin/env python |
---|
2 | |
---|
3 | import unittest, os |
---|
4 | import os.path |
---|
5 | from math import pi, sqrt |
---|
6 | import tempfile |
---|
7 | |
---|
8 | from anuga.config import g, epsilon |
---|
9 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a |
---|
10 | from anuga.utilities.numerical_tools import mean |
---|
11 | from anuga.utilities.polygon import is_inside_polygon |
---|
12 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
13 | from anuga.abstract_2d_finite_volumes.quantity import Quantity |
---|
14 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
15 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
16 | |
---|
17 | from anuga.utilities.system_tools import get_pathname_from_package |
---|
18 | from swb_domain import * |
---|
19 | |
---|
20 | import numpy as num |
---|
21 | |
---|
22 | # Get gateway to C implementation of flux function for direct testing |
---|
23 | from shallow_water_ext import flux_function_central as flux_function |
---|
24 | |
---|
25 | |
---|
26 | |
---|
27 | |
---|
28 | class Test_swb_clean(unittest.TestCase): |
---|
29 | def setUp(self): |
---|
30 | pass |
---|
31 | |
---|
32 | def tearDown(self): |
---|
33 | pass |
---|
34 | |
---|
35 | |
---|
36 | def test_extrema(self): |
---|
37 | """Test that extrema of quantities are computed correctly |
---|
38 | Extrema are updated at every *internal* timestep |
---|
39 | """ |
---|
40 | |
---|
41 | from anuga.abstract_2d_finite_volumes.mesh_factory \ |
---|
42 | import rectangular_cross |
---|
43 | |
---|
44 | initial_runup_height = -0.4 |
---|
45 | final_runup_height = -0.3 |
---|
46 | |
---|
47 | #-------------------------------------------------------------- |
---|
48 | # Setup computational domain |
---|
49 | #-------------------------------------------------------------- |
---|
50 | N = 5 |
---|
51 | points, vertices, boundary = rectangular_cross(N, N) |
---|
52 | domain = Domain(points, vertices, boundary) |
---|
53 | domain.set_name('extrema_test') |
---|
54 | |
---|
55 | #-------------------------------------------------------------- |
---|
56 | # Setup initial conditions |
---|
57 | #-------------------------------------------------------------- |
---|
58 | def topography(x,y): |
---|
59 | return -x/2 # linear bed slope |
---|
60 | |
---|
61 | domain.set_quantity('elevation', topography) # function for elevation |
---|
62 | domain.set_quantity('friction', 0.) # Zero friction |
---|
63 | # Constant negative initial stage |
---|
64 | domain.set_quantity('stage', initial_runup_height) |
---|
65 | domain.set_quantities_to_be_monitored(['stage', 'height'], |
---|
66 | time_interval=[0.5, 2.7], |
---|
67 | polygon=[[0,0], [0,1], |
---|
68 | [1,1], [1,0]]) |
---|
69 | |
---|
70 | assert len(domain.quantities_to_be_monitored) == 2 |
---|
71 | assert domain.quantities_to_be_monitored.has_key('stage') |
---|
72 | assert domain.quantities_to_be_monitored.has_key('height') |
---|
73 | for key in domain.quantities_to_be_monitored['stage'].keys(): |
---|
74 | assert domain.quantities_to_be_monitored['stage'][key] is None |
---|
75 | |
---|
76 | #-------------------------------------------------------------- |
---|
77 | # Setup boundary conditions |
---|
78 | #-------------------------------------------------------------- |
---|
79 | Br = Reflective_boundary(domain) # Reflective wall |
---|
80 | # Constant inflow |
---|
81 | Bd = Dirichlet_boundary([final_runup_height, 0, 0]) |
---|
82 | |
---|
83 | # All reflective to begin with (still water) |
---|
84 | domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) |
---|
85 | |
---|
86 | #-------------------------------------------------------------- |
---|
87 | # Let triangles adjust and check extrema |
---|
88 | #-------------------------------------------------------------- |
---|
89 | for t in domain.evolve(yieldstep=0.1, finaltime=1.0): |
---|
90 | domain.quantity_statistics() # Run it silently |
---|
91 | |
---|
92 | #-------------------------------------------------------------- |
---|
93 | # Test extrema |
---|
94 | #-------------------------------------------------------------- |
---|
95 | stage = domain.quantities_to_be_monitored['stage'] |
---|
96 | assert stage['min'] <= stage['max'] |
---|
97 | |
---|
98 | assert num.allclose(stage['min'], initial_runup_height, |
---|
99 | rtol=1.0/N) # First order accuracy |
---|
100 | |
---|
101 | depth = domain.quantities_to_be_monitored['height'] |
---|
102 | |
---|
103 | assert depth['min'] <= depth['max'] |
---|
104 | assert depth['min'] >= 0.0 |
---|
105 | assert depth['max'] >= 0.0 |
---|
106 | |
---|
107 | #-------------------------------------------------------------- |
---|
108 | # Update boundary to allow inflow |
---|
109 | #-------------------------------------------------------------- |
---|
110 | domain.set_boundary({'right': Bd}) |
---|
111 | |
---|
112 | #-------------------------------------------------------------- |
---|
113 | # Evolve system through time |
---|
114 | #-------------------------------------------------------------- |
---|
115 | for t in domain.evolve(yieldstep=0.1, finaltime=3.0): |
---|
116 | domain.quantity_statistics() # Run it silently |
---|
117 | |
---|
118 | #-------------------------------------------------------------- |
---|
119 | # Test extrema again |
---|
120 | #-------------------------------------------------------------- |
---|
121 | stage = domain.quantities_to_be_monitored['stage'] |
---|
122 | assert stage['min'] <= stage['max'] |
---|
123 | |
---|
124 | assert num.allclose(stage['min'], initial_runup_height, |
---|
125 | rtol = 1.0/N) # First order accuracy |
---|
126 | |
---|
127 | depth = domain.quantities_to_be_monitored['height'] |
---|
128 | assert depth['min'] <= depth['max'] |
---|
129 | assert depth['min'] >= 0.0 |
---|
130 | assert depth['max'] >= 0.0 |
---|
131 | |
---|
132 | # Cleanup |
---|
133 | os.remove(domain.get_name() + '.sww') |
---|
134 | |
---|
135 | |
---|
136 | |
---|
137 | |
---|
138 | |
---|
139 | ################################################################################# |
---|
140 | |
---|
141 | if __name__ == "__main__": |
---|
142 | suite = unittest.makeSuite(Test_swb_clean, 'test') |
---|
143 | runner = unittest.TextTestRunner(verbosity=1) |
---|
144 | runner.run(suite) |
---|