source: anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 7733

Last change on this file since 7733 was 7733, checked in by hudson, 14 years ago

Fixed unit test failures.

  • Property svn:keywords set to LastChangedDate LastChangedRevision LastChangedBy HeadURL Id
File size: 50.3 KB
Line 
1"""
2Finite-volume computations of the shallow water wave equation.
3
4Title: ANGUA shallow_water_domain - 2D triangular domains for finite-volume
5       computations of the shallow water wave equation.
6
7
8Author: Ole Nielsen, Ole.Nielsen@ga.gov.au
9        Stephen Roberts, Stephen.Roberts@anu.edu.au
10        Duncan Gray, Duncan.Gray@ga.gov.au
11
12CreationDate: 2004
13
14Description:
15    This module contains a specialisation of class Domain from
16    module domain.py consisting of methods specific to the
17    Shallow Water Wave Equation
18
19    U_t + E_x + G_y = S
20
21    where
22
23    U = [w, uh, vh]
24    E = [uh, u^2h + gh^2/2, uvh]
25    G = [vh, uvh, v^2h + gh^2/2]
26    S represents source terms forcing the system
27    (e.g. gravity, friction, wind stress, ...)
28
29    and _t, _x, _y denote the derivative with respect to t, x and y
30    respectively.
31
32
33    The quantities are
34
35    symbol    variable name    explanation
36    x         x                horizontal distance from origin [m]
37    y         y                vertical distance from origin [m]
38    z         elevation        elevation of bed on which flow is modelled [m]
39    h         height           water height above z [m]
40    w         stage            absolute water level, w = z+h [m]
41    u                          speed in the x direction [m/s]
42    v                          speed in the y direction [m/s]
43    uh        xmomentum        momentum in the x direction [m^2/s]
44    vh        ymomentum        momentum in the y direction [m^2/s]
45
46    eta                        mannings friction coefficient [to appear]
47    nu                         wind stress coefficient [to appear]
48
49    The conserved quantities are w, uh, vh
50
51Reference:
52    Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
53    Christopher Zoppou and Stephen Roberts,
54    Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
55
56    Hydrodynamic modelling of coastal inundation.
57    Nielsen, O., S. Roberts, D. Gray, A. McPherson and A. Hitchman
58    In Zerger, A. and Argent, R.M. (eds) MODSIM 2005 International Congress on
59    Modelling and Simulation. Modelling and Simulation Society of Australia and
60    New Zealand, December 2005, pp. 518-523. ISBN: 0-9758400-2-9.
61    http://www.mssanz.org.au/modsim05/papers/nielsen.pdf
62
63
64SeeAlso:
65    TRAC administration of ANUGA (User Manuals etc) at
66    https://datamining.anu.edu.au/anuga and Subversion repository at
67    $HeadURL: anuga_core/source/anuga/shallow_water/shallow_water_domain.py $
68
69Constraints: See GPL license in the user guide
70Version: 1.0 ($Revision: 7733 $)
71ModifiedBy:
72    $Author: hudson $
73    $Date: 2010-05-19 00:13:06 +0000 (Wed, 19 May 2010) $
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2010-05-19 00:13:06 +0000 (Wed, 19 May 2010) $
79# $LastChangedRevision: 7733 $
80# $LastChangedBy: hudson $
81
82
83import numpy as num
84
85from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_Domain
86
87from anuga.shallow_water.forcing import Cross_section
88from anuga.pmesh.mesh_interface import create_mesh_from_regions
89from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
90from anuga.geospatial_data.geospatial_data import ensure_geospatial
91
92from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
93
94from anuga.fit_interpolate.interpolate import Modeltime_too_late, \
95                                              Modeltime_too_early
96
97from anuga.geometry.polygon import inside_polygon, polygon_area, \
98                                    is_inside_polygon
99import anuga.utilities.log as log
100
101import types
102from types import IntType, FloatType
103from warnings import warn
104
105
106################################################################################
107# Shallow water domain
108################################################################################
109
110##
111# @brief Class for a shallow water domain.
112class Domain(Generic_Domain):
113
114    #conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
115    #other_quantities = ['elevation', 'friction']
116
117    ##
118    # @brief Instantiate a shallow water domain.
119    # @param coordinates
120    # @param vertices
121    # @param boundary
122    # @param tagged_elements
123    # @param geo_reference
124    # @param use_inscribed_circle
125    # @param mesh_filename
126    # @param use_cache
127    # @param verbose
128    # @param evolved_quantities
129    # @param full_send_dict
130    # @param ghost_recv_dict
131    # @param processor
132    # @param numproc
133    # @param number_of_full_nodes
134    # @param number_of_full_triangles
135    def __init__(self,
136                 coordinates=None,
137                 vertices=None,
138                 boundary=None,
139                 tagged_elements=None,
140                 geo_reference=None,
141                 use_inscribed_circle=False,
142                 mesh_filename=None,
143                 use_cache=False,
144                 verbose=False,
145                 conserved_quantities = None,
146                 evolved_quantities = None,
147                 other_quantities = None,
148                 full_send_dict=None,
149                 ghost_recv_dict=None,
150                 processor=0,
151                 numproc=1,
152                 number_of_full_nodes=None,
153                 number_of_full_triangles=None):
154
155        # Define quantities for the shallow_water domain
156        if conserved_quantities == None:
157            conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
158
159        if evolved_quantities == None:
160            evolved_quantities =  ['stage', 'xmomentum', 'ymomentum']
161           
162        if other_quantities == None:
163            other_quantities = ['elevation', 'friction']
164       
165        Generic_Domain.__init__(self,
166                                coordinates,
167                                vertices,
168                                boundary,
169                                conserved_quantities,
170                                evolved_quantities,
171                                other_quantities,
172                                tagged_elements,
173                                geo_reference,
174                                use_inscribed_circle,
175                                mesh_filename,
176                                use_cache,
177                                verbose,
178                                full_send_dict,
179                                ghost_recv_dict,
180                                processor,
181                                numproc,
182                                number_of_full_nodes=number_of_full_nodes,
183                                number_of_full_triangles=number_of_full_triangles)
184
185        self.set_defaults()
186
187 
188        self.forcing_terms.append(manning_friction_implicit)
189        self.forcing_terms.append(gravity)
190
191        # Stored output
192        self.store = True
193        self.set_store_vertices_uniquely(False)
194
195        self.quantities_to_be_stored = {'elevation': 1, 
196                                        'stage': 2, 
197                                        'xmomentum': 2, 
198                                        'ymomentum': 2}
199
200
201
202
203    ##
204    # @brief Set default values, usually read in from a config file
205    # @param flag
206    def set_defaults(self):
207        """Set the default values in this routine. That way we can inherit class
208        and just over redefine the defaults for the new class
209        """
210
211        from anuga.config import minimum_storable_height
212        from anuga.config import minimum_allowed_height, maximum_allowed_speed
213        from anuga.config import g, epsilon, beta_w, beta_w_dry,\
214             beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
215        from anuga.config import alpha_balance
216        from anuga.config import optimise_dry_cells
217        from anuga.config import optimised_gradient_limiter
218        from anuga.config import use_edge_limiter
219        from anuga.config import use_centroid_velocities
220
221
222
223        self.set_minimum_allowed_height(minimum_allowed_height)
224        self.maximum_allowed_speed = maximum_allowed_speed
225
226        self.g = g
227        self.beta_w = beta_w
228        self.beta_w_dry = beta_w_dry
229        self.beta_uh = beta_uh
230        self.beta_uh_dry = beta_uh_dry
231        self.beta_vh = beta_vh
232        self.beta_vh_dry = beta_vh_dry
233        self.alpha_balance = alpha_balance
234
235        self.tight_slope_limiters = tight_slope_limiters
236        self.optimise_dry_cells = optimise_dry_cells
237
238       
239        self.set_new_mannings_function(False)
240
241        self.minimum_storable_height = minimum_storable_height
242
243         # Limiters
244        self.use_edge_limiter = use_edge_limiter
245        self.optimised_gradient_limiter = optimised_gradient_limiter
246        self.use_centroid_velocities = use_centroid_velocities       
247
248    ##
249    # @brief
250    # @param flag
251    def set_new_mannings_function(self, flag=True):
252        """Cludge to allow unit test to pass, but to
253        also introduce new mannings friction function
254        which takes into account the slope of the bed.
255        The flag is tested in the python wrapper
256        mannings_friction_implicit
257        """
258        if flag:
259            self.use_new_mannings = True
260        else:
261            self.use_new_mannings = False
262
263
264    ##
265    # @brief
266    # @param flag
267    def set_use_edge_limiter(self, flag=True):
268        """Cludge to allow unit test to pass, but to
269        also introduce new edge limiting. The flag is
270        tested in distribute_to_vertices_and_edges
271        """
272        if flag:
273            self.use_edge_limiter = True
274        else:
275            self.use_edge_limiter = False
276
277
278         
279    ##
280    # @brief
281    # @param beta
282    def set_all_limiters(self, beta):
283        """Shorthand to assign one constant value [0,1] to all limiters.
284        0 Corresponds to first order, where as larger values make use of
285        the second order scheme.
286        """
287
288        self.beta_w = beta
289        self.beta_w_dry = beta
290        self.quantities['stage'].beta = beta
291
292        self.beta_uh = beta
293        self.beta_uh_dry = beta
294        self.quantities['xmomentum'].beta = beta
295
296        self.beta_vh = beta
297        self.beta_vh_dry = beta
298        self.quantities['ymomentum'].beta = beta
299
300    ##
301    # @brief
302    # @param flag
303    # @param reduction
304    def set_store_vertices_uniquely(self, flag, reduction=None):
305        """Decide whether vertex values should be stored uniquely as
306        computed in the model (True) or whether they should be reduced to one
307        value per vertex using self.reduction (False).
308        """
309
310        # FIXME (Ole): how about using the word "continuous vertex values" or
311        # "continuous stage surface"
312        self.smooth = not flag
313
314        # Reduction operation for get_vertex_values
315        if reduction is None:
316            self.reduction = mean
317            #self.reduction = min  #Looks better near steep slopes
318
319    ##
320    # @brief Set the minimum depth that will be written to an SWW file.
321    # @param minimum_storable_height The minimum stored height (in m).
322    def set_minimum_storable_height(self, minimum_storable_height):
323        """Set the minimum depth that will be recognised when writing
324        to an sww file. This is useful for removing thin water layers
325        that seems to be caused by friction creep.
326
327        The minimum allowed sww depth is in meters.
328        """
329
330        self.minimum_storable_height = minimum_storable_height
331
332    ##
333    # @brief
334    # @param minimum_allowed_height
335    def set_minimum_allowed_height(self, minimum_allowed_height):
336        """Set minimum depth that will be recognised in the numerical scheme.
337
338        The minimum allowed depth is in meters.
339
340        The parameter H0 (Minimal height for flux computation)
341        is also set by this function
342        """
343
344        #FIXME (Ole): rename H0 to minimum_allowed_height_in_flux_computation
345
346        #FIXME (Ole): Maybe use histogram to identify isolated extreme speeds
347        #and deal with them adaptively similarly to how we used to use 1 order
348        #steps to recover.
349
350        self.minimum_allowed_height = minimum_allowed_height
351        self.H0 = minimum_allowed_height
352
353    ##
354    # @brief
355    # @param maximum_allowed_speed
356    def set_maximum_allowed_speed(self, maximum_allowed_speed):
357        """Set the maximum particle speed that is allowed in water
358        shallower than minimum_allowed_height. This is useful for
359        controlling speeds in very thin layers of water and at the same time
360        allow some movement avoiding pooling of water.
361        """
362
363        self.maximum_allowed_speed = maximum_allowed_speed
364
365    ##
366    # @brief
367    # @param points_file_block_line_size
368    def set_points_file_block_line_size(self, points_file_block_line_size):
369        """Set the minimum depth that will be recognised when writing
370        to an sww file. This is useful for removing thin water layers
371        that seems to be caused by friction creep.
372
373        The minimum allowed sww depth is in meters.
374        """
375        self.points_file_block_line_size = points_file_block_line_size
376
377       
378    # FIXME: Probably obsolete in its curren form   
379    ##
380    # @brief Set the quantities that will be written to an SWW file.
381    # @param q The quantities to be written.
382    # @note Param 'q' may be None, single quantity or list of quantity strings.
383    # @note If 'q' is None, no quantities will be stored in the SWW file.
384    def set_quantities_to_be_stored(self, q):
385        """Specify which quantities will be stored in the sww file
386       
387        q must be either:
388          - a dictionary with quantity names
389          - a list of quantity names (for backwards compatibility)
390          - None
391
392        The format of the dictionary is as follows
393       
394        quantity_name: flag where flag must be either 1 or 2.
395        If flag is 1, the quantity is considered static and will be
396        stored once at the beginning of the simulation in a 1D array.
397       
398        If flag is 2, the quantity is considered time dependent and
399        it will be stored at each yieldstep by appending it to the
400        appropriate 2D array in the sww file.   
401       
402        If q is None, storage will be switched off altogether.
403       
404        Once the simulation has started and thw sww file opened,
405        this function will have no effect.
406       
407        The format, where q is a list of names is for backwards compatibility only.
408        It will take the specified quantities to be time dependent and assume
409        'elevation' to be static regardless.
410        """
411
412        if q is None:
413            self.quantities_to_be_stored = {}
414            self.store = False
415            return
416
417        # Check correcness
418        for quantity_name in q:
419            msg = ('Quantity %s is not a valid conserved quantity'
420                   % quantity_name)
421            assert quantity_name in self.quantities, msg
422
423        if type(q) == types.ListType:
424
425            msg = 'List arguments to set_quantities_to_be_stored '
426            msg += 'has been deprecated and will be removed in future '
427            msg += 'versions of ANUGA.'
428            msg += 'Please use dictionary argument instead'
429            warn(msg, DeprecationWarning) 
430
431       
432       
433            # FIXME: Raise deprecation warning
434            tmp = {}
435            for x in q:
436                tmp[x] = 2
437            tmp['elevation'] = 1   
438            q = tmp     
439           
440        assert type(q) == types.DictType   
441        self.quantities_to_be_stored = q
442
443    ##
444    # @brief
445    # @param indices
446    def get_wet_elements(self, indices=None):
447        """Return indices for elements where h > minimum_allowed_height
448
449        Optional argument:
450            indices is the set of element ids that the operation applies to.
451
452        Usage:
453            indices = get_wet_elements()
454
455        Note, centroid values are used for this operation
456        """
457
458        # Water depth below which it is considered to be 0 in the model
459        # FIXME (Ole): Allow this to be specified as a keyword argument as well
460        from anuga.config import minimum_allowed_height
461
462        elevation = self.get_quantity('elevation').\
463                        get_values(location='centroids', indices=indices)
464        stage = self.get_quantity('stage').\
465                    get_values(location='centroids', indices=indices)
466        depth = stage - elevation
467
468        # Select indices for which depth > 0
469        wet_indices = num.compress(depth > minimum_allowed_height,
470                                   num.arange(len(depth)))
471        return wet_indices
472
473    ##
474    # @brief
475    # @param indices
476    def get_maximum_inundation_elevation(self, indices=None):
477        """Return highest elevation where h > 0
478
479        Optional argument:
480            indices is the set of element ids that the operation applies to.
481
482        Usage:
483            q = get_maximum_inundation_elevation()
484
485        Note, centroid values are used for this operation
486        """
487
488        wet_elements = self.get_wet_elements(indices)
489        return self.get_quantity('elevation').\
490                   get_maximum_value(indices=wet_elements)
491
492    ##
493    # @brief
494    # @param indices
495    def get_maximum_inundation_location(self, indices=None):
496        """Return location of highest elevation where h > 0
497
498        Optional argument:
499            indices is the set of element ids that the operation applies to.
500
501        Usage:
502            q = get_maximum_inundation_location()
503
504        Note, centroid values are used for this operation
505        """
506
507        wet_elements = self.get_wet_elements(indices)
508        return self.get_quantity('elevation').\
509                   get_maximum_location(indices=wet_elements)
510
511
512
513    ##
514    # @brief Get the total flow through an arbitrary poly line.
515    # @param polyline Representation of desired cross section.
516    # @param verbose True if this method is to be verbose.
517    # @note 'polyline' may contain multiple sections allowing complex shapes.
518    # @note Assume absolute UTM coordinates.
519    def get_flow_through_cross_section(self, polyline, verbose=False):
520        """Get the total flow through an arbitrary poly line.
521
522        This is a run-time equivalent of the function with same name
523        in data_manager.py
524
525        Input:
526            polyline: Representation of desired cross section - it may contain
527                      multiple sections allowing for complex shapes. Assume
528                      absolute UTM coordinates.
529                      Format [[x0, y0], [x1, y1], ...]
530
531        Output:
532            Q: Total flow [m^3/s] across given segments.
533        """
534
535
536        cross_section = Cross_section(self, polyline, verbose)
537
538        return cross_section.get_flow_through_cross_section()
539
540
541
542
543    ##
544    # @brief
545    # @param polyline Representation of desired cross section.
546    # @param kind Select energy type to compute ('specific' or 'total').
547    # @param verbose True if this method is to be verbose.
548    # @note 'polyline' may contain multiple sections allowing complex shapes.
549    # @note Assume absolute UTM coordinates.
550    def get_energy_through_cross_section(self, polyline,
551                                         kind='total',
552                                         verbose=False):
553        """Obtain average energy head [m] across specified cross section.
554
555        Inputs:
556            polyline: Representation of desired cross section - it may contain
557                      multiple sections allowing for complex shapes. Assume
558                      absolute UTM coordinates.
559                      Format [[x0, y0], [x1, y1], ...]
560            kind:     Select which energy to compute.
561                      Options are 'specific' and 'total' (default)
562
563        Output:
564            E: Average energy [m] across given segments for all stored times.
565
566        The average velocity is computed for each triangle intersected by
567        the polyline and averaged weighted by segment lengths.
568
569        The typical usage of this function would be to get average energy of
570        flow in a channel, and the polyline would then be a cross section
571        perpendicular to the flow.
572
573        #FIXME (Ole) - need name for this energy reflecting that its dimension
574        is [m].
575        """
576
577
578
579        cross_section = Cross_section(self, polyline, verbose)
580
581        return cross_section.get_energy_through_cross_section(kind)
582
583
584    ##
585    # @brief Run integrity checks on shallow water domain.
586    def check_integrity(self):
587        Generic_Domain.check_integrity(self)
588
589        #Check that we are solving the shallow water wave equation
590        msg = 'First conserved quantity must be "stage"'
591        assert self.conserved_quantities[0] == 'stage', msg
592        msg = 'Second conserved quantity must be "xmomentum"'
593        assert self.conserved_quantities[1] == 'xmomentum', msg
594        msg = 'Third conserved quantity must be "ymomentum"'
595        assert self.conserved_quantities[2] == 'ymomentum', msg
596
597    ##
598    # @brief
599    def extrapolate_second_order_sw(self):
600        #Call correct module function (either from this module or C-extension)
601        extrapolate_second_order_sw(self)
602
603    ##
604    # @brief
605    def compute_fluxes(self):
606        #Call correct module function (either from this module or C-extension)
607        compute_fluxes(self)
608
609    ##
610    # @brief
611    def distribute_to_vertices_and_edges(self):
612        # Call correct module function
613        if self.use_edge_limiter:
614            distribute_using_edge_limiter(self)
615        else:
616            distribute_using_vertex_limiter(self)
617
618
619
620    ##
621    # @brief Evolve the model by one step.
622    # @param yieldstep
623    # @param finaltime
624    # @param duration
625    # @param skip_initial_step
626    def evolve(self,
627               yieldstep=None,
628               finaltime=None,
629               duration=None,
630               skip_initial_step=False):
631        """Specialisation of basic evolve method from parent class"""
632
633        # Call check integrity here rather than from user scripts
634        # self.check_integrity()
635
636        msg = 'Attribute self.beta_w must be in the interval [0, 2]'
637        assert 0 <= self.beta_w <= 2.0, msg
638
639        # Initial update of vertex and edge values before any STORAGE
640        # and or visualisation.
641        # This is done again in the initialisation of the Generic_Domain
642        # evolve loop but we do it here to ensure the values are ok for storage.
643        self.distribute_to_vertices_and_edges()
644
645        if self.store is True and self.time == 0.0:
646            self.initialise_storage()
647
648        # Call basic machinery from parent class
649        for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
650                                       finaltime=finaltime, duration=duration,
651                                       skip_initial_step=skip_initial_step):
652            # Store model data, e.g. for subsequent visualisation
653            if self.store is True:
654                self.store_timestep()
655
656            # Pass control on to outer loop for more specific actions
657            yield(t)
658
659    ##
660    # @brief
661    def initialise_storage(self):
662        """Create and initialise self.writer object for storing data.
663        Also, save x,y and bed elevation
664        """
665
666        from anuga.shallow_water.data_manager import SWW_file
667       
668        # Initialise writer
669        self.writer = SWW_file(self)
670
671        # Store vertices and connectivity
672        self.writer.store_connectivity()
673
674    ##
675    # @brief
676    # @param name
677    def store_timestep(self):
678        """Store time dependent quantities and time.
679
680        Precondition:
681           self.writer has been initialised
682        """
683
684        self.writer.store_timestep()
685
686    ##
687    # @brief Get time stepping statistics string for printing.
688    # @param track_speeds
689    # @param triangle_id
690    def timestepping_statistics(self,
691                                track_speeds=False,
692                                triangle_id=None):
693        """Return string with time stepping statistics for printing or logging
694
695        Optional boolean keyword track_speeds decides whether to report
696        location of smallest timestep as well as a histogram and percentile
697        report.
698        """
699
700        from anuga.config import epsilon, g
701
702        # Call basic machinery from parent class
703        msg = Generic_Domain.timestepping_statistics(self, track_speeds,
704                                                     triangle_id)
705
706        if track_speeds is True:
707            # qwidth determines the text field used for quantities
708            qwidth = self.qwidth
709
710            # Selected triangle
711            k = self.k
712
713            # Report some derived quantities at vertices, edges and centroid
714            # specific to the shallow water wave equation
715            z = self.quantities['elevation']
716            w = self.quantities['stage']
717
718            Vw = w.get_values(location='vertices', indices=[k])[0]
719            Ew = w.get_values(location='edges', indices=[k])[0]
720            Cw = w.get_values(location='centroids', indices=[k])
721
722            Vz = z.get_values(location='vertices', indices=[k])[0]
723            Ez = z.get_values(location='edges', indices=[k])[0]
724            Cz = z.get_values(location='centroids', indices=[k])
725
726            name = 'depth'
727            Vh = Vw-Vz
728            Eh = Ew-Ez
729            Ch = Cw-Cz
730
731            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
732                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
733
734            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
735                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
736
737            s += '    %s: centroid_value = %.4f\n'\
738                 %(name.ljust(qwidth), Ch[0])
739
740            msg += s
741
742            uh = self.quantities['xmomentum']
743            vh = self.quantities['ymomentum']
744
745            Vuh = uh.get_values(location='vertices', indices=[k])[0]
746            Euh = uh.get_values(location='edges', indices=[k])[0]
747            Cuh = uh.get_values(location='centroids', indices=[k])
748
749            Vvh = vh.get_values(location='vertices', indices=[k])[0]
750            Evh = vh.get_values(location='edges', indices=[k])[0]
751            Cvh = vh.get_values(location='centroids', indices=[k])
752
753            # Speeds in each direction
754            Vu = Vuh/(Vh + epsilon)
755            Eu = Euh/(Eh + epsilon)
756            Cu = Cuh/(Ch + epsilon)
757            name = 'U'
758            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
759                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
760
761            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
762                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
763
764            s += '    %s: centroid_value = %.4f\n'\
765                 %(name.ljust(qwidth), Cu[0])
766
767            msg += s
768
769            Vv = Vvh/(Vh + epsilon)
770            Ev = Evh/(Eh + epsilon)
771            Cv = Cvh/(Ch + epsilon)
772            name = 'V'
773            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
774                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
775
776            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
777                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
778
779            s += '    %s: centroid_value = %.4f\n'\
780                 %(name.ljust(qwidth), Cv[0])
781
782            msg += s
783
784            # Froude number in each direction
785            name = 'Froude (x)'
786            Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
787            Efx = Eu/(num.sqrt(g*Eh) + epsilon)
788            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
789
790            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
791                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
792
793            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
794                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
795
796            s += '    %s: centroid_value = %.4f\n'\
797                 %(name.ljust(qwidth), Cfx[0])
798
799            msg += s
800
801            name = 'Froude (y)'
802            Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
803            Efy = Ev/(num.sqrt(g*Eh) + epsilon)
804            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
805
806            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
807                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
808
809            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
810                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
811
812            s += '    %s: centroid_value = %.4f\n'\
813                 %(name.ljust(qwidth), Cfy[0])
814
815            msg += s
816
817        return msg
818       
819       
820
821    def compute_boundary_flows(self):
822        """Compute boundary flows at current timestep.
823                       
824        Quantities computed are:
825           Total inflow across boundary
826           Total outflow across boundary
827           Flow across each tagged boundary segment
828        """
829               
830        # Run through boundary array and compute for each segment
831        # the normal momentum ((uh, vh) dot normal) times segment length.
832        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
833                       
834        # Compute flows along boundary
835       
836        uh = self.get_quantity('xmomentum').get_values(location='edges')
837        vh = self.get_quantity('ymomentum').get_values(location='edges')       
838       
839        # Loop through edges that lie on the boundary and calculate
840        # flows
841        boundary_flows = {}
842        total_boundary_inflow = 0.0
843        total_boundary_outflow = 0.0
844        for vol_id, edge_id in self.boundary:
845            # Compute normal flow across edge. Since normal vector points
846            # away from triangle, a positive sign means that water
847            # flows *out* from this triangle.
848             
849            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
850            normal = self.mesh.get_normal(vol_id, edge_id)
851            length = self.mesh.get_edgelength(vol_id, edge_id)           
852            normal_flow = num.dot(momentum, normal)*length
853           
854            # Reverse sign so that + is taken to mean inflow
855            # and - means outflow. This is more intuitive.
856            edge_flow = -normal_flow
857           
858            # Tally up inflows and outflows separately
859            if edge_flow > 0:
860                # Flow is inflow     
861                total_boundary_inflow += edge_flow                                 
862            else:
863                # Flow is outflow
864                total_boundary_outflow += edge_flow   
865
866            # Tally up flows by boundary tag
867            tag = self.boundary[(vol_id, edge_id)]
868           
869            if tag not in boundary_flows:
870                boundary_flows[tag] = 0.0
871            boundary_flows[tag] += edge_flow
872           
873               
874        return boundary_flows, total_boundary_inflow, total_boundary_outflow
875       
876
877    def compute_forcing_flows(self):
878        """Compute flows in and out of domain due to forcing terms.
879                       
880        Quantities computed are:
881               
882       
883           Total inflow through forcing terms
884           Total outflow through forcing terms
885           Current total volume in domain       
886
887        """
888
889        #FIXME(Ole): We need to separate what part of explicit_update was
890        # due to the normal flux calculations and what is due to forcing terms.
891       
892        pass
893                       
894       
895    def compute_total_volume(self):
896        """Compute total volume (m^3) of water in entire domain
897        """
898       
899        area = self.mesh.get_areas()
900        volume = 0.0
901       
902        stage = self.get_quantity('stage').get_values(location='centroids')
903        elevation = self.get_quantity('elevation').get_values(location='centroids')       
904        depth = stage-elevation
905       
906        return num.sum(depth*area)
907       
908       
909    def volumetric_balance_statistics(self):               
910        """Create volumetric balance report suitable for printing or logging.
911        """
912       
913        (boundary_flows, total_boundary_inflow,
914         total_boundary_outflow) = self.compute_boundary_flows() 
915       
916        s = '---------------------------\n'       
917        s += 'Volumetric balance report:\n'
918        s += '--------------------------\n'
919        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
920        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow       
921        s += 'Net boundary flow by tags [m^3/s]\n'
922        for tag in boundary_flows:
923            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
924       
925        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 
926        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
927       
928        # The go through explicit forcing update and record the rate of change for stage and
929        # record into forcing_inflow and forcing_outflow. Finally compute integral
930        # of depth to obtain total volume of domain.
931       
932        # FIXME(Ole): This part is not yet done.               
933       
934        return s       
935           
936################################################################################
937# End of class Shallow Water Domain
938################################################################################
939
940#-----------------
941# Flux computation
942#-----------------
943
944## @brief Compute fluxes and timestep suitable for all volumes in domain.
945# @param domain The domain to calculate fluxes for.
946def compute_fluxes(domain):
947    """Compute fluxes and timestep suitable for all volumes in domain.
948
949    Compute total flux for each conserved quantity using "flux_function"
950
951    Fluxes across each edge are scaled by edgelengths and summed up
952    Resulting flux is then scaled by area and stored in
953    explicit_update for each of the three conserved quantities
954    stage, xmomentum and ymomentum
955
956    The maximal allowable speed computed by the flux_function for each volume
957    is converted to a timestep that must not be exceeded. The minimum of
958    those is computed as the next overall timestep.
959
960    Post conditions:
961      domain.explicit_update is reset to computed flux values
962      domain.timestep is set to the largest step satisfying all volumes.
963
964    This wrapper calls the underlying C version of compute fluxes
965    """
966
967    import sys
968    from shallow_water_ext import compute_fluxes_ext_central \
969                                  as compute_fluxes_ext
970
971    N = len(domain)    # number_of_triangles
972
973    # Shortcuts
974    Stage = domain.quantities['stage']
975    Xmom = domain.quantities['xmomentum']
976    Ymom = domain.quantities['ymomentum']
977    Bed = domain.quantities['elevation']
978
979    timestep = float(sys.maxint)
980
981    flux_timestep = compute_fluxes_ext(timestep,
982                                       domain.epsilon,
983                                       domain.H0,
984                                       domain.g,
985                                       domain.neighbours,
986                                       domain.neighbour_edges,
987                                       domain.normals,
988                                       domain.edgelengths,
989                                       domain.radii,
990                                       domain.areas,
991                                       domain.tri_full_flag,
992                                       Stage.edge_values,
993                                       Xmom.edge_values,
994                                       Ymom.edge_values,
995                                       Bed.edge_values,
996                                       Stage.boundary_values,
997                                       Xmom.boundary_values,
998                                       Ymom.boundary_values,
999                                       Stage.explicit_update,
1000                                       Xmom.explicit_update,
1001                                       Ymom.explicit_update,
1002                                       domain.already_computed_flux,
1003                                       domain.max_speed,
1004                                       int(domain.optimise_dry_cells))
1005
1006    domain.flux_timestep = flux_timestep
1007
1008################################################################################
1009# Module functions for gradient limiting
1010################################################################################
1011
1012##
1013# @brief Wrapper for C version of extrapolate_second_order_sw.
1014# @param domain The domain to operate on.
1015# @note MH090605 The following method belongs to the shallow_water domain class
1016#       see comments in the corresponding method in shallow_water_ext.c
1017def extrapolate_second_order_sw(domain):
1018    """Wrapper calling C version of extrapolate_second_order_sw"""
1019
1020    import sys
1021    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
1022
1023    N = len(domain) # number_of_triangles
1024
1025    # Shortcuts
1026    Stage = domain.quantities['stage']
1027    Xmom = domain.quantities['xmomentum']
1028    Ymom = domain.quantities['ymomentum']
1029    Elevation = domain.quantities['elevation']
1030
1031    extrapol2(domain,
1032              domain.surrogate_neighbours,
1033              domain.number_of_boundaries,
1034              domain.centroid_coordinates,
1035              Stage.centroid_values,
1036              Xmom.centroid_values,
1037              Ymom.centroid_values,
1038              Elevation.centroid_values,
1039              domain.vertex_coordinates,
1040              Stage.vertex_values,
1041              Xmom.vertex_values,
1042              Ymom.vertex_values,
1043              Elevation.vertex_values,
1044              int(domain.optimise_dry_cells))
1045
1046##
1047# @brief Distribution from centroids to vertices specific to the SWW eqn.
1048# @param domain The domain to operate on.
1049def distribute_using_vertex_limiter(domain):
1050    """Distribution from centroids to vertices specific to the SWW equation.
1051
1052    It will ensure that h (w-z) is always non-negative even in the
1053    presence of steep bed-slopes by taking a weighted average between shallow
1054    and deep cases.
1055
1056    In addition, all conserved quantities get distributed as per either a
1057    constant (order==1) or a piecewise linear function (order==2).
1058
1059    FIXME: more explanation about removal of artificial variability etc
1060
1061    Precondition:
1062      All quantities defined at centroids and bed elevation defined at
1063      vertices.
1064
1065    Postcondition
1066      Conserved quantities defined at vertices
1067    """
1068
1069    # Remove very thin layers of water
1070    protect_against_infinitesimal_and_negative_heights(domain)
1071
1072    # Extrapolate all conserved quantities
1073    if domain.optimised_gradient_limiter:
1074        # MH090605 if second order,
1075        # perform the extrapolation and limiting on
1076        # all of the conserved quantities
1077
1078        if (domain._order_ == 1):
1079            for name in domain.conserved_quantities:
1080                Q = domain.quantities[name]
1081                Q.extrapolate_first_order()
1082        elif domain._order_ == 2:
1083            domain.extrapolate_second_order_sw()
1084        else:
1085            raise 'Unknown order'
1086    else:
1087        # Old code:
1088        for name in domain.conserved_quantities:
1089            Q = domain.quantities[name]
1090
1091            if domain._order_ == 1:
1092                Q.extrapolate_first_order()
1093            elif domain._order_ == 2:
1094                Q.extrapolate_second_order_and_limit_by_vertex()
1095            else:
1096                raise 'Unknown order'
1097
1098    # Take bed elevation into account when water heights are small
1099    balance_deep_and_shallow(domain)
1100
1101    # Compute edge values by interpolation
1102    for name in domain.conserved_quantities:
1103        Q = domain.quantities[name]
1104        Q.interpolate_from_vertices_to_edges()
1105
1106##
1107# @brief Distribution from centroids to edges specific to the SWW eqn.
1108# @param domain The domain to operate on.
1109def distribute_using_edge_limiter(domain):
1110    """Distribution from centroids to edges specific to the SWW eqn.
1111
1112    It will ensure that h (w-z) is always non-negative even in the
1113    presence of steep bed-slopes by taking a weighted average between shallow
1114    and deep cases.
1115
1116    In addition, all conserved quantities get distributed as per either a
1117    constant (order==1) or a piecewise linear function (order==2).
1118
1119
1120    Precondition:
1121      All quantities defined at centroids and bed elevation defined at
1122      vertices.
1123
1124    Postcondition
1125      Conserved quantities defined at vertices
1126    """
1127
1128    # Remove very thin layers of water
1129    protect_against_infinitesimal_and_negative_heights(domain)
1130
1131    for name in domain.conserved_quantities:
1132        Q = domain.quantities[name]
1133        if domain._order_ == 1:
1134            Q.extrapolate_first_order()
1135        elif domain._order_ == 2:
1136            Q.extrapolate_second_order_and_limit_by_edge()
1137        else:
1138            raise 'Unknown order'
1139
1140    balance_deep_and_shallow(domain)
1141
1142    # Compute edge values by interpolation
1143    for name in domain.conserved_quantities:
1144        Q = domain.quantities[name]
1145        Q.interpolate_from_vertices_to_edges()
1146
1147##
1148# @brief  Protect against infinitesimal heights and associated high velocities.
1149# @param domain The domain to operate on.
1150def protect_against_infinitesimal_and_negative_heights(domain):
1151    """Protect against infinitesimal heights and associated high velocities"""
1152
1153    from shallow_water_ext import protect
1154
1155    # Shortcuts
1156    wc = domain.quantities['stage'].centroid_values
1157    zc = domain.quantities['elevation'].centroid_values
1158    xmomc = domain.quantities['xmomentum'].centroid_values
1159    ymomc = domain.quantities['ymomentum'].centroid_values
1160
1161    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1162            domain.epsilon, wc, zc, xmomc, ymomc)
1163
1164##
1165# @brief Wrapper for C function balance_deep_and_shallow_c().
1166# @param domain The domain to operate on.
1167def balance_deep_and_shallow(domain):
1168    """Compute linear combination between stage as computed by
1169    gradient-limiters limiting using w, and stage computed by
1170    gradient-limiters limiting using h (h-limiter).
1171    The former takes precedence when heights are large compared to the
1172    bed slope while the latter takes precedence when heights are
1173    relatively small.  Anything in between is computed as a balanced
1174    linear combination in order to avoid numerical disturbances which
1175    would otherwise appear as a result of hard switching between
1176    modes.
1177
1178    Wrapper for C implementation
1179    """
1180
1181    from shallow_water_ext import balance_deep_and_shallow \
1182                                  as balance_deep_and_shallow_c
1183
1184    # Shortcuts
1185    wc = domain.quantities['stage'].centroid_values
1186    zc = domain.quantities['elevation'].centroid_values
1187    wv = domain.quantities['stage'].vertex_values
1188    zv = domain.quantities['elevation'].vertex_values
1189
1190    # Momentums at centroids
1191    xmomc = domain.quantities['xmomentum'].centroid_values
1192    ymomc = domain.quantities['ymomentum'].centroid_values
1193
1194    # Momentums at vertices
1195    xmomv = domain.quantities['xmomentum'].vertex_values
1196    ymomv = domain.quantities['ymomentum'].vertex_values
1197
1198    balance_deep_and_shallow_c(domain,
1199                               wc, zc, wv, zv, wc,
1200                               xmomc, ymomc, xmomv, ymomv)
1201
1202
1203
1204################################################################################
1205# Standard forcing terms
1206################################################################################
1207
1208##
1209# @brief Apply gravitational pull in the presence of bed slope.
1210# @param domain The domain to apply gravity to.
1211# @note Wrapper for C function gravity_c().
1212def gravity(domain):
1213    """Apply gravitational pull in the presence of bed slope
1214    Wrapper calls underlying C implementation
1215    """
1216
1217    from shallow_water_ext import gravity as gravity_c
1218
1219    xmom_update = domain.quantities['xmomentum'].explicit_update
1220    ymom_update = domain.quantities['ymomentum'].explicit_update
1221
1222    stage = domain.quantities['stage']
1223    elevation = domain.quantities['elevation']
1224
1225    h = stage.centroid_values - elevation.centroid_values
1226    z = elevation.vertex_values
1227
1228    x = domain.get_vertex_coordinates()
1229    g = domain.g
1230
1231    gravity_c(g, h, z, x, xmom_update, ymom_update)    #, 1.0e-6)
1232
1233##
1234# @brief Apply friction to a surface (implicit).
1235# @param domain The domain to apply Manning friction to.
1236# @note Wrapper for C function manning_friction_c().
1237def manning_friction_implicit(domain):
1238    """Apply (Manning) friction to water momentum
1239    Wrapper for c version
1240    """
1241
1242    from shallow_water_ext import manning_friction_old
1243    from shallow_water_ext import manning_friction_new
1244
1245    xmom = domain.quantities['xmomentum']
1246    ymom = domain.quantities['ymomentum']
1247
1248    x = domain.get_vertex_coordinates()
1249   
1250    w = domain.quantities['stage'].centroid_values
1251    z = domain.quantities['elevation'].vertex_values
1252
1253    uh = xmom.centroid_values
1254    vh = ymom.centroid_values
1255    eta = domain.quantities['friction'].centroid_values
1256
1257    xmom_update = xmom.semi_implicit_update
1258    ymom_update = ymom.semi_implicit_update
1259
1260    N = len(domain)
1261    eps = domain.minimum_allowed_height
1262    g = domain.g
1263
1264    if domain.use_new_mannings:
1265        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
1266    else:
1267        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
1268   
1269
1270##
1271# @brief Apply friction to a surface (explicit).
1272# @param domain The domain to apply Manning friction to.
1273# @note Wrapper for C function manning_friction_c().
1274def manning_friction_explicit(domain):
1275    """Apply (Manning) friction to water momentum
1276    Wrapper for c version
1277    """
1278
1279    from shallow_water_ext import manning_friction_old
1280    from shallow_water_ext import manning_friction_new
1281
1282    xmom = domain.quantities['xmomentum']
1283    ymom = domain.quantities['ymomentum']
1284
1285    x = domain.get_vertex_coordinates()
1286   
1287    w = domain.quantities['stage'].centroid_values
1288    z = domain.quantities['elevation'].vertex_values
1289
1290    uh = xmom.centroid_values
1291    vh = ymom.centroid_values
1292    eta = domain.quantities['friction'].centroid_values
1293
1294    xmom_update = xmom.explicit_update
1295    ymom_update = ymom.explicit_update
1296
1297    N = len(domain)
1298    eps = domain.minimum_allowed_height
1299    g = domain.g
1300
1301
1302    if domain.use_new_mannings:
1303        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
1304    else:
1305        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
1306
1307
1308
1309# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1310##
1311# @brief Apply linear friction to a surface.
1312# @param domain The domain to apply Manning friction to.
1313# @note Is this still used (30 Oct 2007)?
1314def linear_friction(domain):
1315    """Apply linear friction to water momentum
1316
1317    Assumes quantity: 'linear_friction' to be present
1318    """
1319
1320    from math import sqrt
1321
1322    w = domain.quantities['stage'].centroid_values
1323    z = domain.quantities['elevation'].centroid_values
1324    h = w-z
1325
1326    uh = domain.quantities['xmomentum'].centroid_values
1327    vh = domain.quantities['ymomentum'].centroid_values
1328    tau = domain.quantities['linear_friction'].centroid_values
1329
1330    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1331    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1332
1333    N = len(domain) # number_of_triangles
1334    eps = domain.minimum_allowed_height
1335    g = domain.g #Not necessary? Why was this added?
1336
1337    for k in range(N):
1338        if tau[k] >= eps:
1339            if h[k] >= eps:
1340                S = -tau[k]/h[k]
1341
1342                #Update momentum
1343                xmom_update[k] += S*uh[k]
1344                ymom_update[k] += S*vh[k]
1345
1346def depth_dependent_friction(domain, default_friction,
1347                             surface_roughness_data,
1348                             verbose=False):
1349    """Returns an array of friction values for each wet element adjusted for depth.
1350
1351    Inputs:
1352        domain - computational domain object
1353        default_friction - depth independent bottom friction
1354        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
1355        friction region.
1356
1357    Outputs:
1358        wet_friction - Array that can be used directly to update friction as follows:
1359                       domain.set_quantity('friction', wet_friction)
1360
1361       
1362       
1363    """
1364
1365    import numpy as num
1366   
1367    # Create a temp array to store updated depth dependent friction for wet elements
1368    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
1369    N=len(domain)
1370    wet_friction    = num.zeros(N, num.float)
1371    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
1372   
1373   
1374    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
1375    # Recompute depth as vector 
1376    d = depth.get_values(location='centroids')
1377 
1378    # rebuild the 'friction' values adjusted for depth at this instant
1379    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
1380       
1381        # Get roughness data for each element
1382        n0 = float(surface_roughness_data[i,0])
1383        d1 = float(surface_roughness_data[i,1])
1384        n1 = float(surface_roughness_data[i,2])
1385        d2 = float(surface_roughness_data[i,3])
1386        n2 = float(surface_roughness_data[i,4])
1387       
1388       
1389        # Recompute friction values from depth for this element
1390               
1391        if d[i]   <= d1: 
1392            depth_dependent_friction = n1
1393        elif d[i] >= d2:
1394            depth_dependent_friction = n2
1395        else:
1396            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
1397           
1398        # check sanity of result
1399        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
1400            log.critical('%s >>>> WARNING: computed depth_dependent friction '
1401                         'out of range, ddf%f, n1=%f, n2=%f'
1402                         % (model_data.basename,
1403                            depth_dependent_friction, n1, n2))
1404       
1405        # update depth dependent friction  for that wet element
1406        wet_friction[i] = depth_dependent_friction
1407       
1408    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
1409   
1410    if verbose :
1411        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
1412        n_min=min(nvals)
1413        n_max=max(nvals)
1414       
1415        log.critical('         ++++ calculate_depth_dependent_friction - '
1416                     'Updated friction - range  %7.3f to %7.3f'
1417                     % (n_min, n_max))
1418   
1419    return wet_friction
1420
1421
1422
1423################################################################################
1424# Initialise module
1425################################################################################
1426
1427from anuga.utilities import compile
1428if compile.can_use_C_extension('shallow_water_ext.c'):
1429    # Underlying C implementations can be accessed
1430    from shallow_water_ext import assign_windfield_values
1431else:
1432    msg = 'C implementations could not be accessed by %s.\n ' % __file__
1433    msg += 'Make sure compile_all.py has been run as described in '
1434    msg += 'the ANUGA installation guide.'
1435    raise Exception, msg
1436
1437# Optimisation with psyco
1438from anuga.config import use_psyco
1439if use_psyco:
1440    try:
1441        import psyco
1442    except:
1443        import os
1444        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
1445            pass
1446            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1447        else:
1448            msg = ('WARNING: psyco (speedup) could not be imported, '
1449                   'you may want to consider installing it')
1450            log.critical(msg)
1451    else:
1452        psyco.bind(Domain.distribute_to_vertices_and_edges)
1453        psyco.bind(Domain.compute_fluxes)
1454
1455
1456if __name__ == "__main__":
1457    pass
Note: See TracBrowser for help on using the repository browser.