source: anuga_core/source/anuga/shallow_water_balanced/test_swb_reporting.py @ 7559

Last change on this file since 7559 was 7559, checked in by steve, 10 years ago

Broke up test_swb_domain so that we can concentrate on new functionality

File size: 5.5 KB
Line 
1#!/usr/bin/env python
2
3import unittest, os
4import os.path
5from math import pi, sqrt
6import tempfile
7
8from anuga.config import g, epsilon
9from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
10from anuga.utilities.numerical_tools import mean
11from anuga.utilities.polygon import is_inside_polygon
12from anuga.coordinate_transforms.geo_reference import Geo_reference
13from anuga.abstract_2d_finite_volumes.quantity import Quantity
14from anuga.geospatial_data.geospatial_data import Geospatial_data
15from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
16
17from anuga.utilities.system_tools import get_pathname_from_package
18from swb_domain import *
19
20import numpy as num
21
22# Get gateway to C implementation of flux function for direct testing
23from shallow_water_ext import flux_function_central as flux_function
24
25
26
27
28class 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', 'stage-elevation'],
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('stage-elevation')
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['stage-elevation']
102        assert depth['min'] <= depth['max']
103        assert depth['min'] >= 0.0
104        assert depth['max'] >= 0.0
105
106        #--------------------------------------------------------------
107        # Update boundary to allow inflow
108        #--------------------------------------------------------------
109        domain.set_boundary({'right': Bd})
110
111        #--------------------------------------------------------------
112        # Evolve system through time
113        #--------------------------------------------------------------
114        for t in domain.evolve(yieldstep=0.1, finaltime=3.0):
115            domain.quantity_statistics()    # Run it silently
116
117        #--------------------------------------------------------------
118        # Test extrema again
119        #--------------------------------------------------------------
120        stage = domain.quantities_to_be_monitored['stage']
121        assert stage['min'] <= stage['max']
122
123        assert num.allclose(stage['min'], initial_runup_height,
124                            rtol = 1.0/N) # First order accuracy
125
126        depth = domain.quantities_to_be_monitored['stage-elevation']
127        assert depth['min'] <= depth['max']
128        assert depth['min'] >= 0.0
129        assert depth['max'] >= 0.0
130
131        # Cleanup
132        os.remove(domain.get_name() + '.sww')
133
134
135
136
137
138#################################################################################
139
140if __name__ == "__main__":
141    suite = unittest.makeSuite(Test_swb_clean, 'test')
142    runner = unittest.TextTestRunner(verbosity=1)
143    runner.run(suite)
Note: See TracBrowser for help on using the repository browser.