source: trunk/anuga_core/source/anuga/shallow_water_balanced/swb_domain.py @ 8072

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

Fixing up swb code

File size: 13.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(2)
83        self.set_default_order(2)
84        self.set_sloped_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        from swb_domain_ext import compute_fluxes_c
122
123        #Shortcuts
124        W  = self.quantities['stage']
125        UH = self.quantities['xmomentum']
126        VH = self.quantities['ymomentum']
127        H  = self.quantities['height']
128        Z  = self.quantities['elevation']
129        U  = self.quantities['xvelocity']
130        V  = self.quantities['yvelocity']
131
132        timestep = self.get_evolve_max_timestep()
133       
134        self.flux_timestep = compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V)
135
136        #print self.flux_timestep
137
138    ##
139    # @brief
140    def distribute_to_vertices_and_edges(self):
141        """Distribution from centroids to edges specific to the SWW eqn.
142
143        It will ensure that h (w-z) is always non-negative even in the
144        presence of steep bed-slopes by taking a weighted average between shallow
145        and deep cases.
146
147        In addition, all conserved quantities get distributed as per either a
148        constant (order==1) or a piecewise linear function (order==2).
149       
150
151        Precondition:
152        All conserved quantities defined at centroids and bed elevation defined at
153        edges.
154       
155        Postcondition
156        Evolved quantities defined at vertices and edges
157        """
158
159        #Sww_domain.distribute_to_vertices_and_edges(self)
160        #return
161   
162        #Shortcuts
163        W  = self.quantities['stage']
164        UH = self.quantities['xmomentum']
165        VH = self.quantities['ymomentum']
166        H  = self.quantities['height']
167        Z  = self.quantities['elevation']
168        U  = self.quantities['xvelocity']
169        V  = self.quantities['yvelocity']
170
171        #Arrays   
172        w_C   = W.centroid_values   
173        uh_C  = UH.centroid_values
174        vh_C  = VH.centroid_values   
175        z_C   = Z.centroid_values
176        h_C   = H.centroid_values
177        u_C   = U.centroid_values
178        v_C   = V.centroid_values
179
180        w_C[:] = num.maximum(w_C, z_C)
181       
182        h_C[:] = w_C - z_C
183
184
185        assert num.min(h_C) >= 0
186               
187        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
188        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
189        num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)       
190
191        H0 = 1.0e-8
192        u_C[:]  = uh_C/(h_C + H0/h_C)
193        v_C[:]  = vh_C/(h_C + H0/h_C)
194
195        num.putmask(h_C, h_C < 1.0e-15, 0.0)
196       
197        for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]:
198            Q = self.quantities[name]
199            if self._order_ == 1:
200                Q.extrapolate_first_order()
201            elif self._order_ == 2:
202                Q.extrapolate_second_order_and_limit_by_edge()
203                #Q.extrapolate_second_order_and_limit_by_vertex()
204            else:
205                raise 'Unknown order'
206
207
208        w_E     = W.edge_values
209        uh_E    = UH.edge_values
210        vh_E    = VH.edge_values       
211        h_E     = H.edge_values
212        z_E     = Z.edge_values
213        u_E     = U.edge_values
214        v_E     = V.edge_values         
215
216
217        #minh_E = num.min(h_E)
218        #msg = 'min h_E = %g ' % minh_E
219        #assert minh_E >= -1.0e-15, msg
220
221        z_E[:]   = w_E - h_E
222
223        num.putmask(h_E, h_E <= 1.0e-15, 0.0)
224        num.putmask(u_E, h_E <= 1.0e-15, 0.0)
225        num.putmask(v_E, h_E <= 1.0e-15, 0.0)
226        num.putmask(w_E, h_E <= 1.0e-15, z_E)
227        #num.putmask(h_E, h_E <= 0.0, 0.0)
228       
229        uh_E[:] = u_E * h_E
230        vh_E[:] = v_E * h_E
231
232        """
233        print '=========================================================='
234        print 'Time ', self.get_time()
235        print h_E
236        print uh_E
237        print vh_E
238        """
239       
240        # Compute vertex values by interpolation
241        for name in self.evolved_quantities:
242            Q = self.quantities[name]
243            Q.interpolate_from_edges_to_vertices()
244
245
246        w_V     = W.vertex_values
247        uh_V    = UH.vertex_values
248        vh_V    = VH.vertex_values     
249        z_V     = Z.vertex_values       
250        h_V     = H.vertex_values
251        u_V     = U.vertex_values
252        v_V     = V.vertex_values               
253
254
255        #w_V[:]    = z_V + h_V
256
257        #num.putmask(u_V, h_V <= 0.0, 0.0)
258        #num.putmask(v_V, h_V <= 0.0, 0.0)
259        #num.putmask(w_V, h_V <= 0.0, z_V)       
260        #num.putmask(h_V, h_V <= 0.0, 0.0)
261       
262        uh_V[:] = u_V * h_V
263        vh_V[:] = v_V * h_V
264
265
266
267
268    ##
269    # @brief
270    def distribute_to_vertices_and_edges_h(self):
271        """Distribution from centroids to edges specific to the SWW eqn.
272
273        It will ensure that h (w-z) is always non-negative even in the
274        presence of steep bed-slopes by taking a weighted average between shallow
275        and deep cases.
276
277        In addition, all conserved quantities get distributed as per either a
278        constant (order==1) or a piecewise linear function (order==2).
279       
280
281        Precondition:
282        All conserved quantities defined at centroids and bed elevation defined at
283        edges.
284       
285        Postcondition
286        Evolved quantities defined at vertices and edges
287        """
288
289
290        #Shortcuts
291        W  = self.quantities['stage']
292        UH = self.quantities['xmomentum']
293        VH = self.quantities['ymomentum']
294        H  = self.quantities['height']
295        Z  = self.quantities['elevation']
296        U  = self.quantities['xvelocity']
297        V  = self.quantities['yvelocity']
298
299        #Arrays   
300        w_C   = W.centroid_values   
301        uh_C  = UH.centroid_values
302        vh_C  = VH.centroid_values   
303        z_C   = Z.centroid_values
304        h_C   = H.centroid_values
305        u_C   = U.centroid_values
306        v_C   = V.centroid_values
307
308        w_C[:] = num.maximum(w_C, z_C)
309       
310        h_C[:] = w_C - z_C
311
312
313        assert num.min(h_C) >= 0
314               
315        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
316        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
317        num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)       
318       
319        u_C[:]  = uh_C/h_C
320        v_C[:]  = vh_C/h_C
321
322        num.putmask(h_C, h_C < 1.0e-15, 0.0)
323       
324        for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]:
325            Q = self.quantities[name]
326            if self._order_ == 1:
327                Q.extrapolate_first_order()
328            elif self._order_ == 2:
329                Q.extrapolate_second_order_and_limit_by_edge()
330                #Q.extrapolate_second_order_and_limit_by_vertex()
331            else:
332                raise 'Unknown order'
333
334
335        w_E     = W.edge_values
336        uh_E    = UH.edge_values
337        vh_E    = VH.edge_values       
338        h_E     = H.edge_values
339        z_E     = Z.edge_values
340        u_E     = U.edge_values
341        v_E     = V.edge_values         
342
343
344        #minh_E = num.min(h_E)
345        #msg = 'min h_E = %g ' % minh_E
346        #assert minh_E >= -1.0e-15, msg
347
348        z_E[:]   = w_E - h_E
349
350        num.putmask(h_E, h_E <= 1.0e-8, 0.0)
351        num.putmask(u_E, h_E <= 1.0e-8, 0.0)
352        num.putmask(v_E, h_E <= 1.0e-8, 0.0)
353        num.putmask(w_E, h_E <= 1.0e-8, z_E)
354        #num.putmask(h_E, h_E <= 0.0, 0.0)
355       
356        uh_E[:] = u_E * h_E
357        vh_E[:] = v_E * h_E
358
359        """
360        print '=========================================================='
361        print 'Time ', self.get_time()
362        print h_E
363        print uh_E
364        print vh_E
365        """
366       
367        # Compute vertex values by interpolation
368        for name in self.evolved_quantities:
369            Q = self.quantities[name]
370            Q.interpolate_from_edges_to_vertices()
371
372
373        w_V     = W.vertex_values
374        uh_V    = UH.vertex_values
375        vh_V    = VH.vertex_values     
376        z_V     = Z.vertex_values       
377        h_V     = H.vertex_values
378        u_V     = U.vertex_values
379        v_V     = V.vertex_values               
380
381
382        #w_V[:]    = z_V + h_V
383
384        #num.putmask(u_V, h_V <= 0.0, 0.0)
385        #num.putmask(v_V, h_V <= 0.0, 0.0)
386        #num.putmask(w_V, h_V <= 0.0, z_V)       
387        #num.putmask(h_V, h_V <= 0.0, 0.0)
388       
389        uh_V[:] = u_V * h_V
390        vh_V[:] = v_V * h_V
391
392
393
394
395    ##
396    # @brief Code to let us use old shallow water domain BCs
397    def conserved_values_to_evolved_values(self, q_cons, q_evol):
398        """Mapping between conserved quantities and the evolved quantities.
399        Used where we have a boundary condition which works with conserved
400        quantities and we now want to use them for the new well balanced
401        code using the evolved quantities
402
403        Typically the old boundary condition will set the values in q_cons,
404
405        q_evol on input will have the values of the evolved quantities at the
406        edge point (useful as it provides values for evlevation).
407        """
408
409        wc  = q_cons[0]
410        uhc = q_cons[1]
411        vhc = q_cons[2]
412
413        we  = q_evol[0]
414        uhe = q_evol[1]
415        vhe = q_evol[2]
416
417        he  = q_evol[3]
418        be  = q_evol[4]
419        ue  = q_evol[5]
420        ve  = q_evol[6]
421
422
423        hc = wc - be
424
425        if hc <= 0.0:
426            hc = 0.0
427            uc = 0.0
428            vc = 0.0
429        else:
430            uc = uhc/hc
431            vc = vhc/hc
432
433        q_evol[0]  = wc
434        q_evol[1]  = uhc
435        q_evol[2]  = vhc
436
437        q_evol[3]  = hc
438        q_evol[4]  = be
439        q_evol[5]  = uc
440        q_evol[6]  = vc
441
442
443        return q_evol
444
445 
Note: See TracBrowser for help on using the repository browser.