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

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

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