source: trunk/anuga_work/shallow_water_balanced_steve/swb_domain.py @ 8248

Last change on this file since 8248 was 8248, checked in by steve, 12 years ago

commiting at the anuga_core level
File size: 14.1 KB
Line 
1
2from anuga.shallow_water.shallow_water_domain import *
3from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain
4
5
6##############################################################################
7# Shallow Water Balanced Domain
8#
9# Uses extra evolved quantities height, elevation, xvelocity, yvelocity
10##############################################################################
11
12class Domain(Sww_domain):
13
14    def __init__(self,
15                 coordinates=None,
16                 vertices=None,
17                 boundary=None,
18                 tagged_elements=None,
19                 geo_reference=None,
20                 use_inscribed_circle=False,
21                 mesh_filename=None,
22                 use_cache=False,
23                 verbose=False,
24                 full_send_dict=None,
25                 ghost_recv_dict=None,
26                 starttime=0.0,
27                 processor=0,
28                 numproc=1,
29                 number_of_full_nodes=None,
30                 number_of_full_triangles=None):
31
32        conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum']
33
34        evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum', \
35                               'height', 'elevation', \
36                               'xvelocity', 'yvelocity' ]
37       
38        other_quantities = [ 'friction', 'x', 'y' ]
39
40       
41        Sww_domain.__init__(self,
42                            coordinates = coordinates,
43                            vertices = vertices,
44                            boundary = boundary,
45                            tagged_elements = tagged_elements,
46                            geo_reference = geo_reference,
47                            use_inscribed_circle = use_inscribed_circle,
48                            mesh_filename = mesh_filename,
49                            use_cache = use_cache,
50                            verbose = verbose,
51                            conserved_quantities = conserved_quantities,
52                            evolved_quantities = evolved_quantities,
53                            other_quantities = other_quantities,
54                            full_send_dict = full_send_dict,
55                            ghost_recv_dict = ghost_recv_dict,
56                            starttime = starttime,
57                            processor = processor,
58                            numproc = numproc,
59                            number_of_full_nodes = number_of_full_nodes,
60                            number_of_full_triangles = number_of_full_triangles)
61       
62        #---------------------
63        # set some defaults
64        #---------------------
65        self.set_timestepping_method(2)
66        self.set_default_order(2)
67        self.set_sloped_mannings_function(True)
68        self.set_centroid_transmissive_bc(True)
69        self.set_CFL(1.0)
70        self.set_beta(1.5)
71        self.quantities['height'].set_beta(1.5)
72
73        #--------------------------------------------
74        # Replace shallow water gravity forcing term
75        # with swb version
76        #--------------------------------------------
77        self.forcing_terms[1] = swb_gravity
78
79        print "Swb_domain"
80
81    def check_integrity(self):
82        Sww_domain.check_integrity(self)
83
84        #Check that the evolved quantities are correct (order)
85        msg = 'First evolved quantity must be "stage"'
86        assert self.evolved_quantities[0] == 'stage', msg
87        msg = 'Second evolved quantity must be "xmomentum"'
88        assert self.evolved_quantities[1] == 'xmomentum', msg
89        msg = 'Third evolved quantity must be "ymomentum"'
90        assert self.evolved_quantities[2] == 'ymomentum', msg
91        msg = 'Fourth evolved quantity must be "height"'
92        assert self.evolved_quantities[3] == 'height', msg
93        msg = 'Fifth evolved quantity must be "elevation"'
94        assert self.evolved_quantities[4] == 'elevation', msg
95        msg = 'Sixth evolved quantity must be "xvelocity"'
96        assert self.evolved_quantities[5] == 'xvelocity', msg       
97        msg = 'Seventh evolved quantity must be "yvelocity"'
98        assert self.evolved_quantities[6] == 'yvelocity', msg       
99
100        msg = 'First other quantity must be "friction"'
101        assert self.other_quantities[0] == 'friction', msg
102        msg = 'Second other quantity must be "x"'
103        assert self.other_quantities[1] == 'x', msg
104        msg = 'Third other quantity must be "y"'
105        assert self.other_quantities[2] == 'y', msg
106
107
108    def compute_fluxes(self):
109        """
110        Call correct module function (either from this module or C-extension)
111        """
112
113        from swb_domain_ext import compute_fluxes_c
114
115        #Shortcuts
116        W  = self.quantities['stage']
117        UH = self.quantities['xmomentum']
118        VH = self.quantities['ymomentum']
119        H  = self.quantities['height']
120        Z  = self.quantities['elevation']
121        U  = self.quantities['xvelocity']
122        V  = self.quantities['yvelocity']
123
124        timestep = self.get_evolve_max_timestep()
125       
126        self.flux_timestep = \
127            compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V)
128
129
130    def distribute_to_vertices_and_edges(self):
131        """Distribution from centroids to edges specific to the SWW eqn.
132
133        It will ensure that h (w-z) is always non-negative even in the
134        presence of steep bed-slopes by taking a weighted average between shallow
135        and deep cases.
136
137        In addition, all conserved quantities get distributed as per either a
138        constant (order==1) or a piecewise linear function (order==2).
139       
140
141        Precondition:
142        All conserved quantities defined at centroids and bed elevation defined at
143        edges.
144       
145        Postcondition
146        Evolved quantities defined at vertices and edges
147        """
148
149        from anuga.shallow_water.shallow_water_domain import \
150                  protect_against_infinitesimal_and_negative_heights as protect
151   
152        #Shortcuts
153        W  = self.quantities['stage']
154        UH = self.quantities['xmomentum']
155        VH = self.quantities['ymomentum']
156        H  = self.quantities['height']
157        Z  = self.quantities['elevation']
158        U  = self.quantities['xvelocity']
159        V  = self.quantities['yvelocity']
160
161        #Arrays   
162        w_C   = W.centroid_values   
163        uh_C  = UH.centroid_values
164        vh_C  = VH.centroid_values   
165        z_C   = Z.centroid_values
166        h_C   = H.centroid_values
167        u_C   = U.centroid_values
168        v_C   = V.centroid_values
169
170        #num_min = num.min(w_C-z_C)
171        #if num_min < -1.0e-5:
172        #    print '**** num.min(w_C-z_C)', num_min
173
174   
175        w_C[:] = num.maximum(w_C, z_C)
176        h_C[:] = w_C - z_C
177
178
179        assert num.min(h_C) >= 0.0
180               
181        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
182        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
183        #num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)
184
185        # Noelle has an alternative method for calculating velcities
186        # Check it out in the GPU shallow water paper
187        H0 = 1.0e-16
188        u_C[:]  = uh_C/(h_C + H0/h_C)
189        v_C[:]  = vh_C/(h_C + H0/h_C)
190
191        #num.putmask(h_C, h_C < 1.0e-15, 0.0)
192       
193        for name in [ 'stage', 'xvelocity', 'yvelocity' ]:
194            Q = self.quantities[name]
195            if self._order_ == 1:
196                Q.extrapolate_first_order()
197            elif self._order_ == 2:
198                Q.extrapolate_second_order_and_limit_by_edge()
199                #Q.extrapolate_second_order_and_limit_by_vertex()
200            else:
201                raise Exception('Unknown order: %s' % str(self._order_))
202
203        for name in [ 'height' ]:
204            Q = self.quantities[name]
205            if self._order_ == 1:
206                Q.extrapolate_first_order()
207            elif self._order_ == 2:
208                #Q.extrapolate_second_order_and_limit_by_edge()
209                Q.extrapolate_second_order_and_limit_by_vertex()
210            else:
211                raise Exception('Unknown order: %s' % str(self._order_))
212
213
214        w_V     = W.vertex_values
215        u_V     = U.vertex_values
216        v_V     = V.vertex_values
217        z_V     = Z.vertex_values
218
219        h_V     = H.vertex_values
220        uh_V    = UH.vertex_values
221        vh_V    = VH.vertex_values
222       
223
224        # Update other quantities
225        #protect(self)
226
227        z_V[:]  = w_V - h_V
228        uh_V[:] = u_V * h_V
229        vh_V[:] = v_V * h_V
230
231       
232        num_min = num.min(h_V)
233        if num_min < -1.0e-14:
234            print 'num.min(h_V)', num_min
235
236       
237        # Compute edge values by interpolation
238        for name in ['xmomentum', 'ymomentum', 'elevation']:
239            Q = self.quantities[name]
240            Q.interpolate_from_vertices_to_edges()
241
242
243
244
245
246    def distribute_to_vertices_and_edges_h(self):
247        """Distribution from centroids to edges specific to the SWW eqn.
248
249        It will ensure that h (w-z) is always non-negative even in the
250        presence of steep bed-slopes by taking a weighted average between shallow
251        and deep cases.
252
253        In addition, all conserved quantities get distributed as per either a
254        constant (order==1) or a piecewise linear function (order==2).
255       
256
257        Precondition:
258        All conserved quantities defined at centroids and bed elevation defined at
259        edges.
260       
261        Postcondition
262        Evolved quantities defined at vertices and edges
263        """
264
265
266        #Shortcuts
267        W  = self.quantities['stage']
268        UH = self.quantities['xmomentum']
269        VH = self.quantities['ymomentum']
270        H  = self.quantities['height']
271        Z  = self.quantities['elevation']
272        U  = self.quantities['xvelocity']
273        V  = self.quantities['yvelocity']
274
275        #Arrays   
276        w_C   = W.centroid_values   
277        uh_C  = UH.centroid_values
278        vh_C  = VH.centroid_values   
279        z_C   = Z.centroid_values
280        h_C   = H.centroid_values
281        u_C   = U.centroid_values
282        v_C   = V.centroid_values
283
284        w_C[:] = num.maximum(w_C, z_C)
285       
286        h_C[:] = w_C - z_C
287
288
289        assert num.min(h_C) >= 0
290               
291        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
292        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
293        num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)       
294       
295        u_C[:]  = uh_C/h_C
296        v_C[:]  = vh_C/h_C
297
298        num.putmask(h_C, h_C < 1.0e-15, 0.0)
299       
300        for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]:
301            Q = self.quantities[name]
302            if self._order_ == 1:
303                Q.extrapolate_first_order()
304            elif self._order_ == 2:
305                Q.extrapolate_second_order_and_limit_by_edge()
306                #Q.extrapolate_second_order_and_limit_by_vertex()
307            else:
308                raise Exception('Unknown order: %s' % str(self._order_))
309
310
311        w_E     = W.edge_values
312        uh_E    = UH.edge_values
313        vh_E    = VH.edge_values       
314        h_E     = H.edge_values
315        z_E     = Z.edge_values
316        u_E     = U.edge_values
317        v_E     = V.edge_values         
318
319
320        #minh_E = num.min(h_E)
321        #msg = 'min h_E = %g ' % minh_E
322        #assert minh_E >= -1.0e-15, msg
323
324        z_E[:]   = w_E - h_E
325
326        num.putmask(h_E, h_E <= 1.0e-8, 0.0)
327        num.putmask(u_E, h_E <= 1.0e-8, 0.0)
328        num.putmask(v_E, h_E <= 1.0e-8, 0.0)
329        num.putmask(w_E, h_E <= 1.0e-8, z_E)
330        #num.putmask(h_E, h_E <= 0.0, 0.0)
331       
332        uh_E[:] = u_E * h_E
333        vh_E[:] = v_E * h_E
334
335        """
336        print '=========================================================='
337        print 'Time ', self.get_time()
338        print h_E
339        print uh_E
340        print vh_E
341        """
342       
343        # Compute vertex values by interpolation
344        for name in self.evolved_quantities:
345            Q = self.quantities[name]
346            Q.interpolate_from_edges_to_vertices()
347
348
349        w_V     = W.vertex_values
350        uh_V    = UH.vertex_values
351        vh_V    = VH.vertex_values     
352        z_V     = Z.vertex_values       
353        h_V     = H.vertex_values
354        u_V     = U.vertex_values
355        v_V     = V.vertex_values               
356
357
358        #w_V[:]    = z_V + h_V
359
360        #num.putmask(u_V, h_V <= 0.0, 0.0)
361        #num.putmask(v_V, h_V <= 0.0, 0.0)
362        #num.putmask(w_V, h_V <= 0.0, z_V)       
363        #num.putmask(h_V, h_V <= 0.0, 0.0)
364       
365        uh_V[:] = u_V * h_V
366        vh_V[:] = v_V * h_V
367
368
369
370
371    def conserved_values_to_evolved_values(self, q_cons, q_evol):
372        """Mapping between conserved quantities and the evolved quantities.
373        Used where we have a boundary condition which works with conserved
374        quantities and we now want to use them for the new well balanced
375        code using the evolved quantities
376
377        Typically the old boundary condition will set the values in q_cons,
378
379        q_evol on input will have the values of the evolved quantities at the
380        edge point (useful as it provides values for evlevation).
381        """
382
383        wc  = q_cons[0]
384        uhc = q_cons[1]
385        vhc = q_cons[2]
386
387        we  = q_evol[0]
388        uhe = q_evol[1]
389        vhe = q_evol[2]
390
391        he  = q_evol[3]
392        be  = q_evol[4]
393        ue  = q_evol[5]
394        ve  = q_evol[6]
395
396
397        hc = wc - be
398
399        if hc <= 0.0:
400            hc = 0.0
401            uc = 0.0
402            vc = 0.0
403        else:
404            uc = uhc/hc
405            vc = vhc/hc
406
407        q_evol[0]  = wc
408        q_evol[1]  = uhc
409        q_evol[2]  = vhc
410
411        q_evol[3]  = hc
412        q_evol[4]  = be
413        q_evol[5]  = uc
414        q_evol[6]  = vc
415
416
417        return q_evol
418
419 
420################################################################################
421# Standard forcing terms
422################################################################################
423
424def swb_gravity(domain):
425    """Apply gravitational pull in the presence of bed slope
426    Wrapper calls underlying C implementation
427    """
428
429    from swb_domain_ext import gravity_c
430
431    #print "Using swb gravity"
432
433    xmom_update = domain.quantities['xmomentum'].explicit_update
434    ymom_update = domain.quantities['ymomentum'].explicit_update
435
436    stage     = domain.quantities['stage']
437    elevation = domain.quantities['elevation']
438
439
440    stage     = stage.vertex_values
441    elevation = elevation.vertex_values
442
443    points = domain.get_vertex_coordinates()
444
445    gravity_c(domain.g, stage, elevation, points, xmom_update, ymom_update)
Note: See TracBrowser for help on using the repository browser.