source: anuga_core/source/anuga/shallow_water_balanced/swb_domain.py @ 7575

Last change on this file since 7575 was 7575, checked in by steve, 13 years ago

Version of compute fluxes with discontinuous bed

File size: 9.5 KB
Line 
1
2from anuga.shallow_water.shallow_water_domain import *
3from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain
4
5from swb_boundary_conditions import Transmissive_boundary
6
7##############################################################################
8# Shallow Water Balanced Domain
9#
10# Uses extra evolved quantities height, elevation, xvelocity, yvelocity
11##############################################################################
12
13##
14# @brief Class for a shallow water balanced domain.
15class Domain(Sww_domain):
16
17    ##
18    # @brief Instantiate a shallow water balanced domain.
19    # @param coordinates
20    # @param vertices
21    # @param boundary
22    # @param tagged_elements
23    # @param geo_reference
24    # @param use_inscribed_circle
25    # @param mesh_filename
26    # @param use_cache
27    # @param verbose
28    # @param full_send_dict
29    # @param ghost_recv_dict
30    # @param processor
31    # @param numproc
32    # @param number_of_full_nodes
33    # @param number_of_full_triangles
34    def __init__(self,
35                 coordinates=None,
36                 vertices=None,
37                 boundary=None,
38                 tagged_elements=None,
39                 geo_reference=None,
40                 use_inscribed_circle=False,
41                 mesh_filename=None,
42                 use_cache=False,
43                 verbose=False,
44                 full_send_dict=None,
45                 ghost_recv_dict=None,
46                 processor=0,
47                 numproc=1,
48                 number_of_full_nodes=None,
49                 number_of_full_triangles=None):
50
51        conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum']
52
53        evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum', \
54                               'height', 'elevation', 'xvelocity', 'yvelocity']
55       
56        other_quantities = [ 'friction' ]
57
58
59        Sww_domain.__init__(self,
60                            coordinates = coordinates,
61                            vertices = vertices,
62                            boundary = boundary,
63                            tagged_elements = tagged_elements,
64                            geo_reference = geo_reference,
65                            use_inscribed_circle = use_inscribed_circle,
66                            mesh_filename = mesh_filename,
67                            use_cache = use_cache,
68                            verbose = verbose,
69                            conserved_quantities = conserved_quantities,
70                            evolved_quantities = evolved_quantities,
71                            other_quantities = other_quantities,
72                            full_send_dict = full_send_dict,
73                            ghost_recv_dict = ghost_recv_dict,
74                            processor = processor,
75                            numproc = numproc,
76                            number_of_full_nodes = number_of_full_nodes,
77                            number_of_full_triangles = number_of_full_triangles)
78       
79        #---------------------
80        # set some defaults
81        #---------------------
82        self.set_timestepping_method(1)
83        self.set_default_order(2)
84        self.set_new_mannings_function(True)
85        self.set_centroid_transmissive_bc(True)
86
87    ##
88    # @brief Run integrity checks on shallow water balanced domain.
89    def check_integrity(self):
90        Sww_domain.check_integrity(self)
91
92        #Check that the evolved quantities are correct (order)
93        msg = 'First evolved quantity must be "stage"'
94        assert self.evolved_quantities[0] == 'stage', msg
95        msg = 'Second evolved quantity must be "xmomentum"'
96        assert self.evolved_quantities[1] == 'xmomentum', msg
97        msg = 'Third evolved quantity must be "ymomentum"'
98        assert self.evolved_quantities[2] == 'ymomentum', msg
99        msg = 'Fourth evolved quantity must be "height"'
100        assert self.evolved_quantities[3] == 'height', msg
101        msg = 'Fifth evolved quantity must be "elevation"'
102        assert self.evolved_quantities[4] == 'elevation', msg
103        msg = 'Sixth evolved quantity must be "xvelocity"'
104        assert self.evolved_quantities[5] == 'xvelocity', msg       
105        msg = 'Seventh evolved quantity must be "yvelocity"'
106        assert self.evolved_quantities[6] == 'yvelocity', msg       
107
108        msg = 'First other quantity must be "friction"'
109        assert self.other_quantities[0] == 'friction', msg
110
111    ##
112    # @brief
113    def compute_fluxes(self):
114        #Call correct module function (either from this module or C-extension)
115
116        from swb_domain_ext import compute_fluxes_c
117
118        #Shortcuts
119        W  = self.quantities['stage']
120        UH = self.quantities['xmomentum']
121        VH = self.quantities['ymomentum']
122        H  = self.quantities['height']
123        Z  = self.quantities['elevation']
124        U  = self.quantities['xvelocity']
125        V  = self.quantities['yvelocity']
126
127        timestep = self.get_evolve_max_timestep()
128       
129        self.flux_timestep = compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V)
130
131    ##
132    # @brief
133    def distribute_to_vertices_and_edges(self):
134        """Distribution from centroids to edges specific to the SWW eqn.
135
136        It will ensure that h (w-z) is always non-negative even in the
137        presence of steep bed-slopes by taking a weighted average between shallow
138        and deep cases.
139
140        In addition, all conserved quantities get distributed as per either a
141        constant (order==1) or a piecewise linear function (order==2).
142       
143
144        Precondition:
145        All conserved quantities defined at centroids and bed elevation defined at
146        edges.
147       
148        Postcondition
149        Evolved quantities defined at vertices and edges
150        """
151
152
153        #Shortcuts
154        W  = self.quantities['stage']
155        UH = self.quantities['xmomentum']
156        VH = self.quantities['ymomentum']
157        H  = self.quantities['height']
158        Z  = self.quantities['elevation']
159        U  = self.quantities['xvelocity']
160        V  = self.quantities['yvelocity']
161
162        #Arrays   
163        w_C   = W.centroid_values   
164        uh_C  = UH.centroid_values
165        vh_C  = VH.centroid_values   
166        z_C   = Z.centroid_values
167        h_C   = H.centroid_values
168        u_C   = U.centroid_values
169        v_C   = V.centroid_values
170
171        w_C[:] = num.maximum(w_C, z_C)
172       
173        h_C[:] = w_C - z_C
174
175
176        assert num.min(h_C) >= 0
177               
178        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
179        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
180        num.putmask(h_C, h_C < 1.0e-15, 1.0e-15)       
181       
182        u_C[:]  = uh_C/h_C
183        v_C[:]  = vh_C/h_C
184       
185        for name in [ 'height', 'elevation', 'xvelocity', 'yvelocity' ]:
186            Q = self.quantities[name]
187            if self._order_ == 1:
188                Q.extrapolate_first_order()
189            elif self._order_ == 2:
190                Q.extrapolate_second_order_and_limit_by_edge()
191            else:
192                raise 'Unknown order'
193
194
195        w_E     = W.edge_values
196        uh_E    = UH.edge_values
197        vh_E    = VH.edge_values       
198        h_E     = H.edge_values
199        z_E     = Z.edge_values
200        u_E     = U.edge_values
201        v_E     = V.edge_values         
202
203
204        w_E[:]   = z_E + h_E
205
206        #num.putmask(u_E, h_temp <= 0.0, 0.0)
207        #num.putmask(v_E, h_temp <= 0.0, 0.0)
208        #num.putmask(w_E, h_temp <= 0.0, z_E+h_E)
209        #num.putmask(h_E, h_E <= 0.0, 0.0)
210       
211        uh_E[:] = u_E * h_E
212        vh_E[:] = v_E * h_E
213
214        """
215        print '=========================================================='
216        print 'Time ', self.get_time()
217        print h_E
218        print uh_E
219        print vh_E
220        """
221       
222        # Compute vertex values by interpolation
223        for name in self.evolved_quantities:
224            Q = self.quantities[name]
225            Q.interpolate_from_edges_to_vertices()
226
227
228        w_V     = W.vertex_values
229        uh_V    = UH.vertex_values
230        vh_V    = VH.vertex_values     
231        z_V     = Z.vertex_values       
232        h_V     = H.vertex_values
233        u_V     = U.vertex_values
234        v_V     = V.vertex_values               
235
236
237        #w_V[:]    = z_V + h_V
238
239        #num.putmask(u_V, h_V <= 0.0, 0.0)
240        #num.putmask(v_V, h_V <= 0.0, 0.0)
241        #num.putmask(w_V, h_V <= 0.0, z_V)       
242        #num.putmask(h_V, h_V <= 0.0, 0.0)
243       
244        uh_V[:] = u_V * h_V
245        vh_V[:] = v_V * h_V
246
247
248    ##
249    # @brief Code to let us use old shallow water domain BCs
250    def conserved_values_to_evolved_values(self, q_cons, q_evol):
251        """Mapping between conserved quantities and the evolved quantities.
252        Used where we have a boundary condition which works with conserved
253        quantities and we now want to use them for the new well balanced
254        code using the evolved quantities
255
256        Typically the old boundary condition will set the values in q_cons,
257
258        q_evol on input will have the values of the evolved quantities at the
259        edge point (useful as it provides values for evlevation).
260        """
261
262        wc  = q_cons[0]
263        uhc = q_cons[1]
264        vhc = q_cons[2]
265
266        we  = q_evol[0]
267        uhe = q_evol[1]
268        vhe = q_evol[2]
269
270        he  = q_evol[3]
271        be  = q_evol[4]
272        ue  = q_evol[5]
273        ve  = q_evol[6]
274
275
276        hc = wc - be
277
278        if hc <= 0.0:
279            hc = 0.0
280            uc = 0.0
281            vc = 0.0
282        else:
283            uc = uhc/hc
284            vc = vhc/hc
285
286        q_evol[0]  = wc
287        q_evol[1]  = uhc
288        q_evol[2]  = vhc
289
290        q_evol[3]  = hc
291        q_evol[4]  = be
292        q_evol[5]  = uc
293        q_evol[6]  = vc
294
295
296        return q_evol
297
298 
Note: See TracBrowser for help on using the repository browser.