source: trunk/anuga_work/development/gareth/balanced_basic/swb2_domain.py @ 8884

Last change on this file since 8884 was 8396, checked in by davies, 13 years ago

balanced_dev: Discouraging use of balanced_basic, and updates to
util.py, and experimental boundary conditions
File size: 14.8 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 -- alternative implementation
8#
9# FIXME: Following the methods in CITE MODSIM PAPER
10#
11##############################################################################
12
13class Domain(Sww_domain):
14
15    def __init__(self,
16                 coordinates=None,
17                 vertices=None,
18                 boundary=None,
19                 tagged_elements=None,
20                 geo_reference=None,
21                 use_inscribed_circle=False,
22                 mesh_filename=None,
23                 use_cache=False,
24                 verbose=False,
25                 full_send_dict=None,
26                 ghost_recv_dict=None,
27                 starttime=0.0,
28                 processor=0,
29                 numproc=1,
30                 number_of_full_nodes=None,
31                 number_of_full_triangles=None):
32
33        conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum']
34
35        evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum']         
36        other_quantities   = [ 'elevation', 'friction', 'height', 
37                               'xvelocity', 'yvelocity', 'x', 'y' ]
38
39       
40        Sww_domain.__init__(self,
41                            coordinates = coordinates,
42                            vertices = vertices,
43                            boundary = boundary,
44                            tagged_elements = tagged_elements,
45                            geo_reference = geo_reference,
46                            use_inscribed_circle = use_inscribed_circle,
47                            mesh_filename = mesh_filename,
48                            use_cache = use_cache,
49                            verbose = verbose,
50                            conserved_quantities = conserved_quantities,
51                            evolved_quantities = evolved_quantities,
52                            other_quantities = other_quantities,
53                            full_send_dict = full_send_dict,
54                            ghost_recv_dict = ghost_recv_dict,
55                            starttime = starttime,
56                            processor = processor,
57                            numproc = numproc,
58                            number_of_full_nodes = number_of_full_nodes,
59                            number_of_full_triangles = number_of_full_triangles)
60       
61        #---------------------
62        # set some defaults
63        #---------------------
64        self.set_CFL(1.0)
65        self.set_use_kinematic_viscosity(False)
66
67        # Because gravity is treated within the flux function,
68        # we remove it from the forcing terms.
69        self.forcing_terms.remove(gravity)
70        print 'Using shallow_water_balanced2 solver in /balanced_basic/'
71        print 'I SUGGEST YOU DO NOT USE THIS ONE -- THE BALANCED DEV SOLVER SEEMS BETTER,'
72        print 'AND IS BETTER TESTED (as of 6/4/2012)'
73        print 'IT PRESENTLY LIVES IN anuga_work/development/gareth/experimental/balanced_dev'
74        print ' '
75        print 'The following error is intended to draw your attention to the message above'
76
77        assert 0==1
78    #-----------------
79    # Flux computation
80    #-----------------
81
82    ## @brief Compute fluxes and timestep suitable for all volumes in domain.
83    # @param domain The domain to calculate fluxes for.
84
85    def compute_fluxes(domain):
86        """Compute fluxes and timestep suitable for all volumes in domain.
87
88        Compute total flux for each conserved quantity using "flux_function"
89
90        Fluxes across each edge are scaled by edgelengths and summed up
91        Resulting flux is then scaled by area and stored in
92        explicit_update for each of the three conserved quantities
93        stage, xmomentum and ymomentum
94
95        The maximal allowable speed computed by the flux_function for each volume
96        is converted to a timestep that must not be exceeded. The minimum of
97        those is computed as the next overall timestep.
98
99        Post conditions:
100          domain.explicit_update is reset to computed flux values
101          domain.timestep is set to the largest step satisfying all volumes.
102
103        This wrapper calls the underlying C version of compute fluxes
104        """
105
106        import sys
107        from swb2_domain_ext import compute_fluxes_ext_central \
108                                      as compute_fluxes_ext
109
110        # Shortcuts
111        Stage = domain.quantities['stage']
112        Xmom = domain.quantities['xmomentum']
113        Ymom = domain.quantities['ymomentum']
114        Bed = domain.quantities['elevation']
115
116        timestep = float(sys.maxint)
117
118        flux_timestep = compute_fluxes_ext(timestep,
119                                           domain.epsilon,
120                                           domain.H0,
121                                           domain.g,
122                                           domain.neighbours,
123                                           domain.neighbour_edges,
124                                           domain.normals,
125                                           domain.edgelengths,
126                                           domain.radii,
127                                           domain.areas,
128                                           domain.tri_full_flag,
129                                           Stage.edge_values,
130                                           Xmom.edge_values,
131                                           Ymom.edge_values,
132                                           Bed.edge_values,
133                                           Stage.boundary_values,
134                                           Xmom.boundary_values,
135                                           Ymom.boundary_values,
136                                           Stage.explicit_update,
137                                           Xmom.explicit_update,
138                                           Ymom.explicit_update,
139                                           domain.already_computed_flux,
140                                           domain.max_speed,
141                                           int(domain.optimise_dry_cells),
142                                           Stage.centroid_values,
143                                           Bed.centroid_values, 
144                                           Bed.vertex_values)
145
146        #import pdb
147        #pdb.set_trace()
148
149        domain.flux_timestep = flux_timestep
150
151
152    def protect_against_infinitesimal_and_negative_heights(domain):
153        """protect against infinitesimal heights and associated high velocities"""
154
155        from swb2_domain_ext import protect
156        #print'using swb2_protect_against ..'
157
158        # shortcuts
159        wc = domain.quantities['stage'].centroid_values
160        zc = domain.quantities['elevation'].centroid_values
161        zv = domain.quantities['elevation'].vertex_values
162        xmomc = domain.quantities['xmomentum'].centroid_values
163        ymomc = domain.quantities['ymomentum'].centroid_values
164
165        protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
166                domain.epsilon, wc, zc,zv, xmomc, ymomc)
167
168    def conserved_values_to_evolved_values(self, q_cons, q_evol):
169        """Mapping between conserved quantities and the evolved quantities.
170        Used where we have a boundary condition which works with conserved
171        quantities and we now want to use them for the new well balanced
172        code using the evolved quantities
173
174        Typically the old boundary condition will set the values in q_cons,
175
176        q_evol on input will have the values of the evolved quantities at the
177        edge point (useful as it provides values for evlevation).
178
179        """
180
181        wc  = q_cons[0]
182        uhc = q_cons[1]
183        vhc = q_cons[2]
184
185
186        q_evol[0]  = wc
187        q_evol[1]  = uhc
188        q_evol[2]  = vhc
189
190        return q_evol
191
192    def distribute_to_vertices_and_edges(self):
193        """ Call correct module function """
194
195                #Shortcuts
196        #W  = self.quantities['stage']
197        #Z  = self.quantities['elevation']
198
199        #Arrays
200        #w_C   = W.centroid_values
201        #z_C   = Z.centroid_values
202
203        #num_min = num.min(w_C-z_C)
204        #if  num_min < 0.0:
205        #    print 'num.min(w_C-z_C)', num_min
206
207        #print 'using swb2 distribute_to_vertices_and_edges'
208        if self.use_edge_limiter:
209            #distribute_using_edge_limiter(self)
210            self.distribute_using_edge_limiter()
211        else:
212            #distribute_using_vertex_limiter(self)
213            self.distribute_using_vertex_limiter()
214
215    #def distribute_to_vertices_and_edges(self):
216    #    """Distribution from centroids to edges specific to the SWW eqn.
217
218    #    In addition, all conserved quantities get distributed as per either a
219    #    constant (order==1) or a piecewise linear function (order==2).
220    #   
221
222    #    Precondition:
223    #    All conserved quantities defined at centroids and bed elevation defined at
224    #    edges.
225    #   
226    #    Postcondition
227    #    Evolved quantities defined at vertices and edges
228    #    """
229    #    from swb2_domain import protect_against_infinitesimal_and_negative_heights
230
231    #    if self._order_ == 1:
232    #        for name in [ 'stage', 'xmomentum', 'ymomentum' ]:
233    #            Q = self.quantities[name]
234    #            Q.extrapolate_first_order()
235    #    elif self._order_ == 2:
236    #        for name in [ 'stage', 'xmomentum', 'ymomentum' ]:
237    #            Q = self.quantities[name]
238    #            #Q.extrapolate_second_order_and_limit_by_edge()
239    #            Q.extrapolate_second_order_and_limit_by_vertex()
240
241    #        #self.extrapolate_second_order_sw()
242    #    else:
243    #        raise Exception('Unknown order: %s' % str(self._order_))
244
245
246    #    # Update other quantities
247    #    protect_against_infinitesimal_and_negative_heights(self)
248
249    #   
250    #    # Compute edge values by interpolation
251    #    for name in ['stage', 'xmomentum', 'ymomentum']:
252    #        Q = self.quantities[name]
253    #        Q.interpolate_from_vertices_to_edges()
254
255
256    def distribute_using_vertex_limiter(domain):
257        """Distribution from centroids to vertices specific to the SWW equation.
258
259        Precondition:
260          All quantities defined at centroids and bed elevation defined at
261          vertices.
262
263        Postcondition
264          Conserved quantities defined at vertices
265        """
266        #from swb2_domain import protect_against_infinitesimal_and_negative_heights
267
268        #print 'using swb2 distribute_using_vertex_limiter'
269        # Remove very thin layers of water
270        domain.protect_against_infinitesimal_and_negative_heights()
271
272        # Extrapolate all conserved quantities
273        if domain.optimised_gradient_limiter:
274            # MH090605 if second order,
275            # perform the extrapolation and limiting on
276            # all of the conserved quantities
277
278            if (domain._order_ == 1):
279                for name in domain.conserved_quantities:
280                    Q = domain.quantities[name]
281                    Q.extrapolate_first_order()
282            elif domain._order_ == 2:
283                domain.extrapolate_second_order_sw()
284            else:
285                raise Exception('Unknown order')
286        else:
287            # Old code:
288            for name in domain.conserved_quantities:
289                Q = domain.quantities[name]
290
291                if domain._order_ == 1:
292                    Q.extrapolate_first_order()
293                elif domain._order_ == 2:
294                    Q.extrapolate_second_order_and_limit_by_vertex()
295                else:
296                    raise Exception('Unknown order')
297
298        # Take bed elevation into account when water heights are small
299        #balance_deep_and_shallow(domain)
300
301        # Compute edge values by interpolation
302        for name in domain.conserved_quantities:
303            Q = domain.quantities[name]
304            Q.interpolate_from_vertices_to_edges()
305
306
307    def distribute_using_edge_limiter(domain):
308        """Distribution from centroids to edges specific to the SWW eqn.
309
310        Precondition:
311          All quantities defined at centroids and bed elevation defined at
312          vertices.
313
314        Postcondition
315          Conserved quantities defined at vertices
316        """
317        #from swb2_domain import protect_against_infinitesimal_and_negative_heights
318
319        # Remove very thin layers of water
320        domain.protect_against_infinitesimal_and_negative_heights()
321
322        #print 'using swb2 distribute_using_vertex_limiter'
323        for name in domain.conserved_quantities:
324            Q = domain.quantities[name]
325            if domain._order_ == 1:
326                Q.extrapolate_first_order()
327            elif domain._order_ == 2:
328                Q.extrapolate_second_order_and_limit_by_edge()
329            else:
330                raise Exception('Unknown order')
331
332        #balance_deep_and_shallow(domain)
333
334        # Compute edge values by interpolation
335        for name in domain.conserved_quantities:
336            Q = domain.quantities[name]
337            Q.interpolate_from_vertices_to_edges()
338
339
340    def update_centroids_of_velocities_and_height(self):
341        """Calculate the centroid values of velocities and height based
342        on the values of the quantities stage and x and y momentum
343
344        Assumes that stage and momentum are up to date
345
346        Useful for kinematic viscosity calculations
347        """
348
349        # For shallow water we need to update height xvelocity and yvelocity
350
351        #Shortcuts
352        W  = self.quantities['stage']
353        UH = self.quantities['xmomentum']
354        VH = self.quantities['ymomentum']
355        H  = self.quantities['height']
356        Z  = self.quantities['elevation']
357        U  = self.quantities['xvelocity']
358        V  = self.quantities['yvelocity']
359
360        #print num.min(W.centroid_values)
361
362        # Make sure boundary values of conserved quantites
363        # are consistent with value of functions at centroids
364        #self.distribute_to_vertices_and_edges()
365        Z.set_boundary_values_from_edges()
366
367        #W.set_boundary_values_from_edges()
368        #UH.set_boundary_values_from_edges()
369        #VH.set_boundary_values_from_edges()
370
371        # Update height values
372        H.set_values(W.centroid_values-Z.centroid_values, location='centroids')
373        H.set_boundary_values( num.where(W.boundary_values-Z.boundary_values>=0,
374                                         W.boundary_values-Z.boundary_values, 0.0))
375
376        #assert num.min(H.centroid_values) >= 0
377        #assert num.min(H.boundary_values) >= 0
378
379        #Aliases
380        uh_C  = UH.centroid_values
381        vh_C  = VH.centroid_values
382        h_C   = H.centroid_values
383
384        uh_B  = UH.boundary_values
385        vh_B  = VH.boundary_values
386        h_B   = H.boundary_values
387
388        H0 = 1.0e-8
389       
390        U.set_values(uh_C/(h_C + H0/h_C), location='centroids')
391        V.set_values(vh_C/(h_C + H0/h_C), location='centroids')
392
393        U.set_boundary_values(uh_B/(h_B + H0/h_B))
394        V.set_boundary_values(vh_B/(h_B + H0/h_B))
Note: See TracBrowser for help on using the repository browser.