source: anuga_work/development/sudi/sw_2d/swb_domain.py @ 7741

Last change on this file since 7741 was 7741, checked in by mungkasi, 15 years ago

Changing to be consistent with old anuga_1_1 code
File size: 10.0 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(2)
83        self.set_default_order(2)
84        self.set_new_mannings_function(True)
85        self.set_centroid_transmissive_bc(True)
86        self.set_CFL(1.0)
87        self.set_beta(1.0)
88
89    ##
90    # @brief Run integrity checks on shallow water balanced domain.
91    def check_integrity(self):
92        Sww_domain.check_integrity(self)
93
94        #Check that the evolved quantities are correct (order)
95        msg = 'First evolved quantity must be "stage"'
96        assert self.evolved_quantities[0] == 'stage', msg
97        msg = 'Second evolved quantity must be "xmomentum"'
98        assert self.evolved_quantities[1] == 'xmomentum', msg
99        msg = 'Third evolved quantity must be "ymomentum"'
100        assert self.evolved_quantities[2] == 'ymomentum', msg
101        msg = 'Fourth evolved quantity must be "height"'
102        assert self.evolved_quantities[3] == 'height', msg
103        msg = 'Fifth evolved quantity must be "elevation"'
104        assert self.evolved_quantities[4] == 'elevation', msg
105        msg = 'Sixth evolved quantity must be "xvelocity"'
106        assert self.evolved_quantities[5] == 'xvelocity', msg       
107        msg = 'Seventh evolved quantity must be "yvelocity"'
108        assert self.evolved_quantities[6] == 'yvelocity', msg       
109
110        msg = 'First other quantity must be "friction"'
111        assert self.other_quantities[0] == 'friction', msg
112
113    ##
114    # @brief
115    def compute_fluxes(self):
116        #Call correct module function (either from this module or C-extension)
117
118        compute_fluxes(self)
119
120        return
121
122
123   
124        from swb_domain_ext import compute_fluxes_c
125
126        #Shortcuts
127        W  = self.quantities['stage']
128        UH = self.quantities['xmomentum']
129        VH = self.quantities['ymomentum']
130        H  = self.quantities['height']
131        Z  = self.quantities['elevation']
132        U  = self.quantities['xvelocity']
133        V  = self.quantities['yvelocity']
134
135        timestep = self.get_evolve_max_timestep()
136       
137        self.stable_timestep = compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V)
138
139        #print self.stable_timestep
140
141    ##
142    # @brief
143    def distribute_to_vertices_and_edges(self):
144        """Distribution from centroids to edges specific to the SWW eqn.
145
146        It will ensure that h (w-z) is always non-negative even in the
147        presence of steep bed-slopes by taking a weighted average between shallow
148        and deep cases.
149
150        In addition, all conserved quantities get distributed as per either a
151        constant (order==1) or a piecewise linear function (order==2).
152       
153
154        Precondition:
155        All conserved quantities defined at centroids and bed elevation defined at
156        edges.
157       
158        Postcondition
159        Evolved quantities defined at vertices and edges
160        """
161
162        Sww_domain.distribute_to_vertices_and_edges(self)
163        return
164
165   
166        #Shortcuts
167        W  = self.quantities['stage']
168        UH = self.quantities['xmomentum']
169        VH = self.quantities['ymomentum']
170        H  = self.quantities['height']
171        Z  = self.quantities['elevation']
172        U  = self.quantities['xvelocity']
173        V  = self.quantities['yvelocity']
174
175        #Arrays   
176        w_C   = W.centroid_values   
177        uh_C  = UH.centroid_values
178        vh_C  = VH.centroid_values   
179        z_C   = Z.centroid_values
180        h_C   = H.centroid_values
181        u_C   = U.centroid_values
182        v_C   = V.centroid_values
183
184        w_C[:] = num.maximum(w_C, z_C)
185       
186        h_C[:] = w_C - z_C
187
188
189        assert num.min(h_C) >= 0
190               
191        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
192        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
193        num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)       
194
195        H0 = 1.0e-16
196        u_C[:]  = uh_C/(h_C + H0/h_C)
197        v_C[:]  = vh_C/(h_C + H0/h_C)
198
199        num.putmask(h_C, h_C < 1.0e-15, 0.0)
200       
201        for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]:
202            Q = self.quantities[name]
203            if self._order_ == 1:
204                Q.extrapolate_first_order()
205            elif self._order_ == 2:
206                Q.extrapolate_second_order_and_limit_by_edge()
207                #Q.extrapolate_second_order_and_limit_by_vertex()
208            else:
209                raise 'Unknown order'
210
211
212        w_E     = W.edge_values
213        uh_E    = UH.edge_values
214        vh_E    = VH.edge_values       
215        h_E     = H.edge_values
216        z_E     = Z.edge_values
217        u_E     = U.edge_values
218        v_E     = V.edge_values         
219
220
221        #minh_E = num.min(h_E)
222        #msg = 'min h_E = %g ' % minh_E
223        #assert minh_E >= -1.0e-15, msg
224
225        z_E[:]   = w_E - h_E
226
227        num.putmask(h_E, h_E <= 1.0e-15, 0.0)
228        num.putmask(u_E, h_E <= 1.0e-15, 0.0)
229        num.putmask(v_E, h_E <= 1.0e-15, 0.0)
230        num.putmask(w_E, h_E <= 1.0e-15, z_E)
231        #num.putmask(h_E, h_E <= 0.0, 0.0)
232       
233        uh_E[:] = u_E * h_E
234        vh_E[:] = v_E * h_E
235
236        """
237        print '=========================================================='
238        print 'Time ', self.get_time()
239        print h_E
240        print uh_E
241        print vh_E
242        """
243       
244        # Compute vertex values by interpolation
245        for name in self.evolved_quantities:
246            Q = self.quantities[name]
247            Q.interpolate_from_edges_to_vertices()
248
249
250        w_V     = W.vertex_values
251        uh_V    = UH.vertex_values
252        vh_V    = VH.vertex_values     
253        z_V     = Z.vertex_values       
254        h_V     = H.vertex_values
255        u_V     = U.vertex_values
256        v_V     = V.vertex_values               
257
258
259        #w_V[:]    = z_V + h_V
260
261        #num.putmask(u_V, h_V <= 0.0, 0.0)
262        #num.putmask(v_V, h_V <= 0.0, 0.0)
263        #num.putmask(w_V, h_V <= 0.0, z_V)       
264        #num.putmask(h_V, h_V <= 0.0, 0.0)
265       
266        uh_V[:] = u_V * h_V
267        vh_V[:] = v_V * h_V
268
269
270
271
272
273    ##
274    # @brief Code to let us use old shallow water domain BCs
275    def conserved_values_to_evolved_values(self, q_cons, q_evol):
276        """Mapping between conserved quantities and the evolved quantities.
277        Used where we have a boundary condition which works with conserved
278        quantities and we now want to use them for the new well balanced
279        code using the evolved quantities
280
281        Typically the old boundary condition will set the values in q_cons,
282
283        q_evol on input will have the values of the evolved quantities at the
284        edge point (useful as it provides values for evlevation).
285        """
286
287        wc  = q_cons[0]
288        uhc = q_cons[1]
289        vhc = q_cons[2]
290
291        we  = q_evol[0]
292        uhe = q_evol[1]
293        vhe = q_evol[2]
294
295        he  = q_evol[3]
296        be  = q_evol[4]
297        ue  = q_evol[5]
298        ve  = q_evol[6]
299
300
301        hc = wc - be
302
303        if hc <= 0.0:
304            hc = 0.0
305            uc = 0.0
306            vc = 0.0
307        else:
308            uc = uhc/hc
309            vc = vhc/hc
310
311        q_evol[0]  = wc
312        q_evol[1]  = uhc
313        q_evol[2]  = vhc
314
315        q_evol[3]  = hc
316        q_evol[4]  = be
317        q_evol[5]  = uc
318        q_evol[6]  = vc
319
320
321        return q_evol
322
323 
Note: See TracBrowser for help on using the repository browser.