source: trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py @ 8356

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

Balanced_dev minor efficiency improvements
File size: 20.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 -- 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        # Most of these override the options in config.py
64        #------------------------------------------------
65        self.set_CFL(1.0)
66        self.set_use_kinematic_viscosity(False)
67        self.timestepping_method='rk2' 
68        # The following allows storage of the negative depths associated with this method
69        self.minimum_storable_height=-99999999999.0 
70
71        self.use_edge_limiter=True
72        self.default_order=2
73        self.extrapolate_velocity_second_order=True
74
75        # Note that the extrapolation method used in quantity_ext.c (e.g.
76        # extrapolate_second_order_and_limit_by_edge) uses a constant value for
77        # all the betas == beta_rk2 in our case. 
78        self.beta_w=2.0
79        self.beta_w_dry=2.0
80        self.beta_uh=2.0
81        self.beta_uh_dry=2.0
82        self.beta_vh=2.0
83        self.beta_vh_dry=2.0
84
85        #self.optimise_dry_cells=True
86        # "self.optimise_dry_cells=False" presently ensures that the stage is
87        # always >= minimum_bed_edge value.  Actually, we should only need to
88        # apply 'False' on the very first time step (to deal with stage values
89        # that were initialised below the bed by the user). After this, the
90        # algorithm should take care of itself, and 'True' should be okay.
91        self.optimise_dry_cells=False 
92
93        # Because gravity is treated within the flux function,
94        # we remove it from the forcing terms.
95        self.forcing_terms.remove(gravity)
96
97        # We need the edge_coordinates for the extrapolation
98        self.edge_coordinates=self.get_edge_midpoint_coordinates()
99
100        # We demand that vertex values are stored uniquely
101        self.set_store_vertices_smoothly(False)
102
103        print '##########################################################################'
104        print '#'
105        print '# Using shallow_water_balanced2 solver in experimental/balanced_dev/'
106        print "#"
107        print "# This solver is experimental - it may cause your computer to explode "
108        print "#"
109        print "# It is not expected to perform well in problems with very"
110        print "# shallow water flowing down steep slopes (such that the stage_centroid_value "
111        print "# is less than the maximum bed_edge_value on a given triangle) "
112        print "#"
113        print "# Also, it allows the stage_centroid_value to drop to the minimum bed_edge_value"
114        print "# on it's triangle. In other ANUGA versions (e.g. 1.2.1), the limit would be the"
115        print "# bed_centroid_value. This means that triangles store slightly more water than they are"
116        print "# typically interpreted to store, which might have significance in some applications."
117        print "#"
118        print "# You will probably be able to tell if the above are causing you problems by convergence testing"
119        print "#"
120        print '# Note that many options in config.py have been overridden by the solver'
121        print '#'
122        print '##########################################################################'
123
124    #-------------------------------
125    # Specialisation of Evolve
126    #-------------------------------
127    # This alters the 'evolve' method from shallow_water_domain.py so that just
128    # before the file-output step, the vertex values are replaced with the
129    # centroid values.
130    # This is useful so that we can unambigously know what the centroid values
131    # of various parameters are
132
133    def evolve(self,
134               yieldstep=None,
135               finaltime=None,
136               duration=None,
137               skip_initial_step=False):
138        """Specialisation of basic evolve method from parent class.
139       
140            Evolve the model by 1 step.
141        """
142
143        # Call check integrity here rather than from user scripts
144        # self.check_integrity()
145
146        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
147        assert 0 <= self.beta_w <= 2.0, msg
148
149        # Initial update of vertex and edge values before any STORAGE
150        # and or visualisation.
151        # This is done again in the initialisation of the Generic_Domain
152        # evolve loop but we do it here to ensure the values are ok for storage.
153        self.distribute_to_vertices_and_edges()
154
155        if self.store is True and self.get_time() == self.get_starttime():
156            self.initialise_storage()
157
158        # Call basic machinery from parent class
159        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
160                                       finaltime=finaltime, duration=duration,
161                                       skip_initial_step=skip_initial_step):
162            # Store model data, e.g. for subsequent visualisation
163            if self.store is True:
164                # FIXME: Extrapolation is done before writing, because I had
165                # trouble correctly computing the centroid values from the
166                # vertex values (an 'average' did not seem to work correctly in
167                # very shallow cells -- there was a discrepency between
168                # domain.quantity['blah'].centroid_values and the values
169                # computed from the sww using the vertex averge). There should
170                # be a much much more disk-efficient way to store the centroid
171                # values than this
172                self.extrapolate_second_order_edge_sw()
173
174                # Store the timestep
175                self.store_timestep()
176
177            # Pass control on to outer loop for more specific actions
178            yield(t)
179   
180    #-----------------
181    # Flux computation
182    #-----------------
183
184    ## @brief Compute fluxes and timestep suitable for all volumes in domain.
185    # @param domain The domain to calculate fluxes for.
186    def compute_fluxes(domain):
187        """Compute fluxes and timestep suitable for all volumes in domain.
188
189        Compute total flux for each conserved quantity using "flux_function"
190
191        Fluxes across each edge are scaled by edgelengths and summed up
192        Resulting flux is then scaled by area and stored in
193        explicit_update for each of the three conserved quantities
194        stage, xmomentum and ymomentum
195
196        The maximal allowable speed computed by the flux_function for each volume
197        is converted to a timestep that must not be exceeded. The minimum of
198        those is computed as the next overall timestep.
199
200        Post conditions:
201          domain.explicit_update is reset to computed flux values
202          domain.timestep is set to the largest step satisfying all volumes.
203
204        This wrapper calls the underlying C version of compute fluxes
205        """
206
207        import sys
208        from swb2_domain_ext import compute_fluxes_ext_central \
209                                      as compute_fluxes_ext
210
211        # Shortcuts
212        Stage = domain.quantities['stage']
213        Xmom = domain.quantities['xmomentum']
214        Ymom = domain.quantities['ymomentum']
215        Bed = domain.quantities['elevation']
216
217        timestep = float(sys.maxint)
218
219        flux_timestep = compute_fluxes_ext(timestep,
220                                           domain.epsilon,
221                                           domain.H0,
222                                           domain.g,
223                                           domain.neighbours,
224                                           domain.neighbour_edges,
225                                           domain.normals,
226                                           domain.edgelengths,
227                                           domain.radii,
228                                           domain.areas,
229                                           domain.tri_full_flag,
230                                           Stage.edge_values,
231                                           Xmom.edge_values,
232                                           Ymom.edge_values,
233                                           Bed.edge_values,
234                                           Stage.boundary_values,
235                                           Xmom.boundary_values,
236                                           Ymom.boundary_values,
237                                           Stage.explicit_update,
238                                           Xmom.explicit_update,
239                                           Ymom.explicit_update,
240                                           domain.already_computed_flux,
241                                           domain.max_speed,
242                                           int(domain.optimise_dry_cells),
243                                           Stage.centroid_values,
244                                           Bed.centroid_values, 
245                                           Bed.vertex_values)
246
247        #import pdb
248        #pdb.set_trace()
249
250        domain.flux_timestep = flux_timestep
251
252
253    def protect_against_infinitesimal_and_negative_heights(domain):
254        """protect against infinitesimal heights and associated high velocities"""
255
256        from swb2_domain_ext import protect
257        #print'using swb2_protect_against ..'
258
259        # shortcuts
260        wc = domain.quantities['stage'].centroid_values
261        wv = domain.quantities['stage'].vertex_values
262        zc = domain.quantities['elevation'].centroid_values
263        zv = domain.quantities['elevation'].vertex_values
264        xmomc = domain.quantities['xmomentum'].centroid_values
265        ymomc = domain.quantities['ymomentum'].centroid_values
266
267        protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
268                domain.epsilon, wc, wv, zc,zv, xmomc, ymomc)
269
270    def conserved_values_to_evolved_values(self, q_cons, q_evol):
271        """Mapping between conserved quantities and the evolved quantities.
272        Used where we have a boundary condition which works with conserved
273        quantities and we now want to use them for the new well balanced
274        code using the evolved quantities
275
276        Typically the old boundary condition will set the values in q_cons,
277
278        q_evol on input will have the values of the evolved quantities at the
279        edge point (useful as it provides values for evlevation).
280
281        """
282
283        wc  = q_cons[0]
284        uhc = q_cons[1]
285        vhc = q_cons[2]
286
287
288        q_evol[0]  = wc
289        q_evol[1]  = uhc
290        q_evol[2]  = vhc
291
292        return q_evol
293
294    def distribute_to_vertices_and_edges(self):
295        """ Call correct module function """
296
297        if self.use_edge_limiter:
298            #distribute_using_edge_limiter(self)
299            self.distribute_using_edge_limiter()
300        else:
301            #distribute_using_vertex_limiter(self)
302            self.distribute_using_vertex_limiter()
303
304
305    def distribute_using_vertex_limiter(domain):
306        """Distribution from centroids to vertices specific to the SWW equation.
307
308        Precondition:
309          All quantities defined at centroids and bed elevation defined at
310          vertices.
311
312        Postcondition
313          Conserved quantities defined at vertices
314        """
315        # Remove very thin layers of water
316        domain.protect_against_infinitesimal_and_negative_heights()
317
318        # Extrapolate all conserved quantities
319        if domain.optimised_gradient_limiter:
320            # MH090605 if second order,
321            # perform the extrapolation and limiting on
322            # all of the conserved quantities
323
324            if (domain._order_ == 1):
325                for name in domain.conserved_quantities:
326                    Q = domain.quantities[name]
327                    Q.extrapolate_first_order()
328            elif domain._order_ == 2:
329                domain.extrapolate_second_order_sw()
330            else:
331                raise Exception('Unknown order')
332        else:
333            # Old code:
334            for name in domain.conserved_quantities:
335                Q = domain.quantities[name]
336
337                if domain._order_ == 1:
338                    Q.extrapolate_first_order()
339                elif domain._order_ == 2:
340                    Q.extrapolate_second_order_and_limit_by_vertex()
341                else:
342                    raise Exception('Unknown order')
343
344        # Take bed elevation into account when water heights are small
345        #balance_deep_and_shallow(domain)
346
347        # Compute edge values by interpolation
348        for name in domain.conserved_quantities:
349            Q = domain.quantities[name]
350            Q.interpolate_from_vertices_to_edges()
351
352
353    def distribute_using_edge_limiter(domain):
354        """Distribution from centroids to edges specific to the SWW eqn.
355
356        Precondition:
357          All quantities defined at centroids and bed elevation defined at
358          vertices.
359
360        Postcondition
361          Conserved quantities defined at vertices
362        """
363        #from swb2_domain import protect_against_infinitesimal_and_negative_heights
364
365        # Remove very thin layers of water
366        domain.protect_against_infinitesimal_and_negative_heights()
367
368        #for name in domain.conserved_quantities:
369        #    Q = domain.quantities[name]
370        #    if domain._order_ == 1:
371        #        Q.extrapolate_first_order()
372        #    elif domain._order_ == 2:
373        #        #Q.extrapolate_second_order_and_limit_by_edge()
374        #        # FIXME: This use of 'break' is hacky
375        #        domain.extrapolate_second_order_edge_sw()
376        #        break
377        #    else:
378        #        raise Exception('Unknown order')
379       
380        domain.extrapolate_second_order_edge_sw()
381
382        #balance_deep_and_shallow(domain)
383
384        # Compute vertex values by interpolation
385        #for name in domain.conserved_quantities:
386        #    Q = domain.quantities[name]
387        #    Q.interpolate_from_edges_to_vertices()
388        #    #Q.interpolate_from_vertices_to_edges()
389
390
391    def update_other_quantities(self):
392        """ There may be a need to calculates some of the other quantities
393        based on the new values of conserved quantities
394
395        But that is not needed in this version of the solver.
396        """
397        pass
398        # The centroid values of height and x and y velocity
399        # might not have been setup
400       
401        #self.update_centroids_of_velocities_and_height()
402        #
403        #for name in ['height', 'xvelocity', 'yvelocity']:
404        #    Q = self.quantities[name]
405        #    if self._order_ == 1:
406        #        Q.extrapolate_first_order()
407        #    elif self._order_ == 2:
408        #        Q.extrapolate_second_order_and_limit_by_edge()
409        #    else:
410        #        raise Exception('Unknown order')
411
412
413    def update_centroids_of_velocities_and_height(self):
414        """Calculate the centroid values of velocities and height based
415        on the values of the quantities stage and x and y momentum
416
417        Assumes that stage and momentum are up to date
418
419        Useful for kinematic viscosity calculations
420        """
421        pass
422        ## For shallow water we need to update height xvelocity and yvelocity
423
424        ##Shortcuts
425        #W  = self.quantities['stage']
426        #UH = self.quantities['xmomentum']
427        #VH = self.quantities['ymomentum']
428        #H  = self.quantities['height']
429        #Z  = self.quantities['elevation']
430        #U  = self.quantities['xvelocity']
431        #V  = self.quantities['yvelocity']
432
433        ##print num.min(W.centroid_values)
434
435        ## Make sure boundary values of conserved quantites
436        ## are consistent with value of functions at centroids
437        ##self.distribute_to_vertices_and_edges()
438        #Z.set_boundary_values_from_edges()
439
440        ##W.set_boundary_values_from_edges()
441        ##UH.set_boundary_values_from_edges()
442        ##VH.set_boundary_values_from_edges()
443
444        ## Update height values
445        #H.set_values(W.centroid_values-Z.centroid_values, location='centroids')
446        #H.set_boundary_values( num.where(W.boundary_values-Z.boundary_values>=0,
447        #                                 W.boundary_values-Z.boundary_values, 0.0))
448
449        ##assert num.min(H.centroid_values) >= 0
450        ##assert num.min(H.boundary_values) >= 0
451
452        ##Aliases
453        #uh_C  = UH.centroid_values
454        #vh_C  = VH.centroid_values
455        #h_C   = H.centroid_values
456
457        #uh_B  = UH.boundary_values
458        #vh_B  = VH.boundary_values
459        #h_B   = H.boundary_values
460
461        #H0 = 1.0e-8
462        #
463        #U.set_values(uh_C/(h_C + H0/h_C), location='centroids')
464        #V.set_values(vh_C/(h_C + H0/h_C), location='centroids')
465
466        #U.set_boundary_values(uh_B/(h_B + H0/h_B))
467        #V.set_boundary_values(vh_B/(h_B + H0/h_B))
468
469
470
471    def extrapolate_second_order_edge_sw(self):
472        """Wrapper calling C version of extrapolate_second_order_sw, using
473            edges instead of vertices in the extrapolation step. The routine
474            then interpolates from edges to vertices.
475            The momentum extrapolation / interpolation is based on either
476            momentum or velocity, depending on the choice of
477            extrapolate_velocity_second_order.
478
479        self  the domain to operate on
480
481        """
482
483        from swb2_domain_ext import extrapolate_second_order_edge_sw as extrapol2
484
485        # Shortcuts
486        Stage = self.quantities['stage']
487        Xmom = self.quantities['xmomentum']
488        Ymom = self.quantities['ymomentum']
489        Elevation = self.quantities['elevation']
490
491        extrapol2(self,
492                  self.surrogate_neighbours,
493                  self.number_of_boundaries,
494                  self.centroid_coordinates,
495                  Stage.centroid_values,
496                  Xmom.centroid_values,
497                  Ymom.centroid_values,
498                  Elevation.centroid_values,
499                  self.edge_coordinates,
500                  Stage.edge_values,
501                  Xmom.edge_values,
502                  Ymom.edge_values,
503                  Elevation.edge_values,
504                  Stage.vertex_values,
505                  Xmom.vertex_values,
506                  Ymom.vertex_values,
507                  Elevation.vertex_values,
508                  int(self.optimise_dry_cells),
509                  int(self.extrapolate_velocity_second_order))
Note: See TracBrowser for help on using the repository browser.