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

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

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