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

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

Still working on channel flow
File size: 17.3 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_w_v(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
206
207        w_V     = W.vertex_values
208        u_V     = U.vertex_values
209        v_V     = V.vertex_values
210        z_V     = Z.vertex_values
211
212        h_V     = H.vertex_values
213        uh_V    = UH.vertex_values
214        vh_V    = VH.vertex_values
215       
216
217        # Update other quantities
218        #protect(self)
219
220        h_V[:]  = w_V - z_V
221        uh_V[:] = u_V * h_V
222        vh_V[:] = v_V * h_V
223
224       
225        num_min = num.min(h_V)
226        if num_min < -1.0e-14:
227            #print 'num.min(h_V)', num_min
228            pass
229
230       
231        # Compute edge values by interpolation
232        for name in ['xmomentum', 'ymomentum', 'elevation']:
233            Q = self.quantities[name]
234            Q.interpolate_from_vertices_to_edges()
235
236
237
238
239    def distribute_to_vertices_and_edges_w_h(self):
240        """Distribution from centroids to edges specific to the SWW eqn.
241
242        It will ensure that h (w-z) is always non-negative even in the
243        presence of steep bed-slopes by taking a weighted average between shallow
244        and deep cases.
245
246        In addition, all conserved quantities get distributed as per either a
247        constant (order==1) or a piecewise linear function (order==2).
248
249
250        Precondition:
251        All conserved quantities defined at centroids and bed elevation defined at
252        edges.
253
254        Postcondition
255        Evolved quantities defined at vertices and edges
256        """
257
258        from anuga.shallow_water.shallow_water_domain import \
259                  protect_against_infinitesimal_and_negative_heights as protect
260
261        #Shortcuts
262        W  = self.quantities['stage']
263        UH = self.quantities['xmomentum']
264        VH = self.quantities['ymomentum']
265        H  = self.quantities['height']
266        Z  = self.quantities['elevation']
267        U  = self.quantities['xvelocity']
268        V  = self.quantities['yvelocity']
269
270        #Arrays
271        w_C   = W.centroid_values
272        uh_C  = UH.centroid_values
273        vh_C  = VH.centroid_values
274        z_C   = Z.centroid_values
275        h_C   = H.centroid_values
276        u_C   = U.centroid_values
277        v_C   = V.centroid_values
278
279        #num_min = num.min(w_C-z_C)
280        #if num_min < -1.0e-5:
281        #    print '**** num.min(w_C-z_C)', num_min
282
283
284        w_C[:] = num.maximum(w_C, z_C)
285        h_C[:] = w_C - z_C
286
287
288        assert num.min(h_C) >= 0.0
289
290        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
291        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
292        #num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)
293
294        # Noelle has an alternative method for calculating velcities
295        # Check it out in the GPU shallow water paper
296        H0 = 1.0e-16
297        u_C[:]  = uh_C/(h_C + H0/h_C)
298        v_C[:]  = vh_C/(h_C + H0/h_C)
299
300        #num.putmask(h_C, h_C < 1.0e-15, 0.0)
301
302        for name in [ 'stage', 'xvelocity', 'yvelocity' ]:
303            Q = self.quantities[name]
304            if self._order_ == 1:
305                Q.extrapolate_first_order()
306            elif self._order_ == 2:
307                Q.extrapolate_second_order_and_limit_by_edge()
308                #Q.extrapolate_second_order_and_limit_by_vertex()
309            else:
310                raise Exception('Unknown order: %s' % str(self._order_))
311
312        for name in [ 'height' ]:
313            Q = self.quantities[name]
314            if self._order_ == 1:
315                Q.extrapolate_first_order()
316            elif self._order_ == 2:
317                #Q.extrapolate_second_order_and_limit_by_edge()
318                Q.extrapolate_second_order_and_limit_by_vertex()
319            else:
320                raise Exception('Unknown order: %s' % str(self._order_))
321
322
323        w_V     = W.vertex_values
324        u_V     = U.vertex_values
325        v_V     = V.vertex_values
326        z_V     = Z.vertex_values
327
328        h_V     = H.vertex_values
329        uh_V    = UH.vertex_values
330        vh_V    = VH.vertex_values
331
332
333        # Update other quantities
334        #protect(self)
335
336        z_V[:]  = w_V - h_V
337        uh_V[:] = u_V * h_V
338        vh_V[:] = v_V * h_V
339
340
341        num_min = num.min(h_V)
342        if num_min < -1.0e-14:
343            #print 'num.min(h_V)', num_min
344            pass
345
346
347        # Compute edge values by interpolation
348        for name in ['xmomentum', 'ymomentum', 'elevation']:
349            Q = self.quantities[name]
350            Q.interpolate_from_vertices_to_edges()
351
352
353
354
355    def distribute_to_vertices_and_edges_h(self):
356        """Distribution from centroids to edges specific to the SWW eqn.
357
358        It will ensure that h (w-z) is always non-negative even in the
359        presence of steep bed-slopes by taking a weighted average between shallow
360        and deep cases.
361
362        In addition, all conserved quantities get distributed as per either a
363        constant (order==1) or a piecewise linear function (order==2).
364       
365
366        Precondition:
367        All conserved quantities defined at centroids and bed elevation defined at
368        edges.
369       
370        Postcondition
371        Evolved quantities defined at vertices and edges
372        """
373
374
375        #Shortcuts
376        W  = self.quantities['stage']
377        UH = self.quantities['xmomentum']
378        VH = self.quantities['ymomentum']
379        H  = self.quantities['height']
380        Z  = self.quantities['elevation']
381        U  = self.quantities['xvelocity']
382        V  = self.quantities['yvelocity']
383
384        #Arrays   
385        w_C   = W.centroid_values   
386        uh_C  = UH.centroid_values
387        vh_C  = VH.centroid_values   
388        z_C   = Z.centroid_values
389        h_C   = H.centroid_values
390        u_C   = U.centroid_values
391        v_C   = V.centroid_values
392
393        w_C[:] = num.maximum(w_C, z_C)
394       
395        h_C[:] = w_C - z_C
396
397
398        assert num.min(h_C) >= 0
399               
400        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
401        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
402        num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)       
403       
404        u_C[:]  = uh_C/h_C
405        v_C[:]  = vh_C/h_C
406
407        num.putmask(h_C, h_C < 1.0e-15, 0.0)
408       
409        for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]:
410            Q = self.quantities[name]
411            if self._order_ == 1:
412                Q.extrapolate_first_order()
413            elif self._order_ == 2:
414                Q.extrapolate_second_order_and_limit_by_edge()
415                #Q.extrapolate_second_order_and_limit_by_vertex()
416            else:
417                raise Exception('Unknown order: %s' % str(self._order_))
418
419
420        w_E     = W.edge_values
421        uh_E    = UH.edge_values
422        vh_E    = VH.edge_values       
423        h_E     = H.edge_values
424        z_E     = Z.edge_values
425        u_E     = U.edge_values
426        v_E     = V.edge_values         
427
428
429        #minh_E = num.min(h_E)
430        #msg = 'min h_E = %g ' % minh_E
431        #assert minh_E >= -1.0e-15, msg
432
433        z_E[:]   = w_E - h_E
434
435        num.putmask(h_E, h_E <= 1.0e-8, 0.0)
436        num.putmask(u_E, h_E <= 1.0e-8, 0.0)
437        num.putmask(v_E, h_E <= 1.0e-8, 0.0)
438        num.putmask(w_E, h_E <= 1.0e-8, z_E)
439        #num.putmask(h_E, h_E <= 0.0, 0.0)
440       
441        uh_E[:] = u_E * h_E
442        vh_E[:] = v_E * h_E
443
444        """
445        print '=========================================================='
446        print 'Time ', self.get_time()
447        print h_E
448        print uh_E
449        print vh_E
450        """
451       
452        # Compute vertex values by interpolation
453        for name in self.evolved_quantities:
454            Q = self.quantities[name]
455            Q.interpolate_from_edges_to_vertices()
456
457
458        w_V     = W.vertex_values
459        uh_V    = UH.vertex_values
460        vh_V    = VH.vertex_values     
461        z_V     = Z.vertex_values       
462        h_V     = H.vertex_values
463        u_V     = U.vertex_values
464        v_V     = V.vertex_values               
465
466
467        #w_V[:]    = z_V + h_V
468
469        #num.putmask(u_V, h_V <= 0.0, 0.0)
470        #num.putmask(v_V, h_V <= 0.0, 0.0)
471        #num.putmask(w_V, h_V <= 0.0, z_V)       
472        #num.putmask(h_V, h_V <= 0.0, 0.0)
473       
474        uh_V[:] = u_V * h_V
475        vh_V[:] = v_V * h_V
476
477
478
479
480    def conserved_values_to_evolved_values(self, q_cons, q_evol):
481        """Mapping between conserved quantities and the evolved quantities.
482        Used where we have a boundary condition which works with conserved
483        quantities and we now want to use them for the new well balanced
484        code using the evolved quantities
485
486        Typically the old boundary condition will set the values in q_cons,
487
488        q_evol on input will have the values of the evolved quantities at the
489        edge point (useful as it provides values for evlevation).
490        """
491
492        wc  = q_cons[0]
493        uhc = q_cons[1]
494        vhc = q_cons[2]
495
496        we  = q_evol[0]
497        uhe = q_evol[1]
498        vhe = q_evol[2]
499
500        he  = q_evol[3]
501        be  = q_evol[4]
502        ue  = q_evol[5]
503        ve  = q_evol[6]
504
505
506        hc = wc - be
507
508        if hc <= 0.0:
509            hc = 0.0
510            uc = 0.0
511            vc = 0.0
512        else:
513            uc = uhc/hc
514            vc = vhc/hc
515
516        q_evol[0]  = wc
517        q_evol[1]  = uhc
518        q_evol[2]  = vhc
519
520        q_evol[3]  = hc
521        q_evol[4]  = be
522        q_evol[5]  = uc
523        q_evol[6]  = vc
524
525
526        return q_evol
527
528 
529################################################################################
530# Standard forcing terms
531################################################################################
532
533def swb_gravity(domain):
534    """Apply gravitational pull in the presence of bed slope
535    Wrapper calls underlying C implementation
536    """
537
538    from swb_domain_ext import gravity_c
539
540    #print "Using swb gravity"
541
542    xmom_update = domain.quantities['xmomentum'].explicit_update
543    ymom_update = domain.quantities['ymomentum'].explicit_update
544
545    stage     = domain.quantities['stage']
546    elevation = domain.quantities['elevation']
547
548
549    stage     = stage.vertex_values
550    elevation = elevation.vertex_values
551
552    points = domain.get_vertex_coordinates()
553
554    gravity_c(domain.g, stage, elevation, points, xmom_update, ymom_update)
Note: See TracBrowser for help on using the repository browser.