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

Last change on this file since 8302 was 8302, checked in by davies, 12 years ago

Adding files for swb2

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