source: trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py @ 7776

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

Removed redundant data_manager class. Unit tests are running, but may fail.

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