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

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

Shallow water refactorings - all unit tests pass, and new file conversion test module created.

  • 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: 7736 $)
71ModifiedBy:
72    $Author: hudson $
73    $Date: 2010-05-24 02:30:34 +0000 (Mon, 24 May 2010) $
74"""
75
76# Subversion keywords:
77#
78# $LastChangedDate: 2010-05-24 02:30:34 +0000 (Mon, 24 May 2010) $
79# $LastChangedRevision: 7736 $
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
92
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        from anuga.shallow_water.sww_file import SWW_file
668       
669        # Initialise writer
670        self.writer = SWW_file(self)
671
672        # Store vertices and connectivity
673        self.writer.store_connectivity()
674
675    ##
676    # @brief
677    # @param name
678    def store_timestep(self):
679        """Store time dependent quantities and time.
680
681        Precondition:
682           self.writer has been initialised
683        """
684
685        self.writer.store_timestep()
686
687    ##
688    # @brief Get time stepping statistics string for printing.
689    # @param track_speeds
690    # @param triangle_id
691    def timestepping_statistics(self,
692                                track_speeds=False,
693                                triangle_id=None):
694        """Return string with time stepping statistics for printing or logging
695
696        Optional boolean keyword track_speeds decides whether to report
697        location of smallest timestep as well as a histogram and percentile
698        report.
699        """
700
701        from anuga.config import epsilon, g
702
703        # Call basic machinery from parent class
704        msg = Generic_Domain.timestepping_statistics(self, track_speeds,
705                                                     triangle_id)
706
707        if track_speeds is True:
708            # qwidth determines the text field used for quantities
709            qwidth = self.qwidth
710
711            # Selected triangle
712            k = self.k
713
714            # Report some derived quantities at vertices, edges and centroid
715            # specific to the shallow water wave equation
716            z = self.quantities['elevation']
717            w = self.quantities['stage']
718
719            Vw = w.get_values(location='vertices', indices=[k])[0]
720            Ew = w.get_values(location='edges', indices=[k])[0]
721            Cw = w.get_values(location='centroids', indices=[k])
722
723            Vz = z.get_values(location='vertices', indices=[k])[0]
724            Ez = z.get_values(location='edges', indices=[k])[0]
725            Cz = z.get_values(location='centroids', indices=[k])
726
727            name = 'depth'
728            Vh = Vw-Vz
729            Eh = Ew-Ez
730            Ch = Cw-Cz
731
732            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
733                 %(name.ljust(qwidth), Vh[0], Vh[1], Vh[2])
734
735            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
736                 %(name.ljust(qwidth), Eh[0], Eh[1], Eh[2])
737
738            s += '    %s: centroid_value = %.4f\n'\
739                 %(name.ljust(qwidth), Ch[0])
740
741            msg += s
742
743            uh = self.quantities['xmomentum']
744            vh = self.quantities['ymomentum']
745
746            Vuh = uh.get_values(location='vertices', indices=[k])[0]
747            Euh = uh.get_values(location='edges', indices=[k])[0]
748            Cuh = uh.get_values(location='centroids', indices=[k])
749
750            Vvh = vh.get_values(location='vertices', indices=[k])[0]
751            Evh = vh.get_values(location='edges', indices=[k])[0]
752            Cvh = vh.get_values(location='centroids', indices=[k])
753
754            # Speeds in each direction
755            Vu = Vuh/(Vh + epsilon)
756            Eu = Euh/(Eh + epsilon)
757            Cu = Cuh/(Ch + epsilon)
758            name = 'U'
759            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
760                 %(name.ljust(qwidth), Vu[0], Vu[1], Vu[2])
761
762            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
763                 %(name.ljust(qwidth), Eu[0], Eu[1], Eu[2])
764
765            s += '    %s: centroid_value = %.4f\n'\
766                 %(name.ljust(qwidth), Cu[0])
767
768            msg += s
769
770            Vv = Vvh/(Vh + epsilon)
771            Ev = Evh/(Eh + epsilon)
772            Cv = Cvh/(Ch + epsilon)
773            name = 'V'
774            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
775                 %(name.ljust(qwidth), Vv[0], Vv[1], Vv[2])
776
777            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
778                 %(name.ljust(qwidth), Ev[0], Ev[1], Ev[2])
779
780            s += '    %s: centroid_value = %.4f\n'\
781                 %(name.ljust(qwidth), Cv[0])
782
783            msg += s
784
785            # Froude number in each direction
786            name = 'Froude (x)'
787            Vfx = Vu/(num.sqrt(g*Vh) + epsilon)
788            Efx = Eu/(num.sqrt(g*Eh) + epsilon)
789            Cfx = Cu/(num.sqrt(g*Ch) + epsilon)
790
791            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
792                 %(name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])
793
794            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
795                 %(name.ljust(qwidth), Efx[0], Efx[1], Efx[2])
796
797            s += '    %s: centroid_value = %.4f\n'\
798                 %(name.ljust(qwidth), Cfx[0])
799
800            msg += s
801
802            name = 'Froude (y)'
803            Vfy = Vv/(num.sqrt(g*Vh) + epsilon)
804            Efy = Ev/(num.sqrt(g*Eh) + epsilon)
805            Cfy = Cv/(num.sqrt(g*Ch) + epsilon)
806
807            s  = '    %s: vertex_values =  %.4f,\t %.4f,\t %.4f\n'\
808                 %(name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])
809
810            s += '    %s: edge_values =    %.4f,\t %.4f,\t %.4f\n'\
811                 %(name.ljust(qwidth), Efy[0], Efy[1], Efy[2])
812
813            s += '    %s: centroid_value = %.4f\n'\
814                 %(name.ljust(qwidth), Cfy[0])
815
816            msg += s
817
818        return msg
819       
820       
821
822    def compute_boundary_flows(self):
823        """Compute boundary flows at current timestep.
824                       
825        Quantities computed are:
826           Total inflow across boundary
827           Total outflow across boundary
828           Flow across each tagged boundary segment
829        """
830               
831        # Run through boundary array and compute for each segment
832        # the normal momentum ((uh, vh) dot normal) times segment length.
833        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
834                       
835        # Compute flows along boundary
836       
837        uh = self.get_quantity('xmomentum').get_values(location='edges')
838        vh = self.get_quantity('ymomentum').get_values(location='edges')       
839       
840        # Loop through edges that lie on the boundary and calculate
841        # flows
842        boundary_flows = {}
843        total_boundary_inflow = 0.0
844        total_boundary_outflow = 0.0
845        for vol_id, edge_id in self.boundary:
846            # Compute normal flow across edge. Since normal vector points
847            # away from triangle, a positive sign means that water
848            # flows *out* from this triangle.
849             
850            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
851            normal = self.mesh.get_normal(vol_id, edge_id)
852            length = self.mesh.get_edgelength(vol_id, edge_id)           
853            normal_flow = num.dot(momentum, normal)*length
854           
855            # Reverse sign so that + is taken to mean inflow
856            # and - means outflow. This is more intuitive.
857            edge_flow = -normal_flow
858           
859            # Tally up inflows and outflows separately
860            if edge_flow > 0:
861                # Flow is inflow     
862                total_boundary_inflow += edge_flow                                 
863            else:
864                # Flow is outflow
865                total_boundary_outflow += edge_flow   
866
867            # Tally up flows by boundary tag
868            tag = self.boundary[(vol_id, edge_id)]
869           
870            if tag not in boundary_flows:
871                boundary_flows[tag] = 0.0
872            boundary_flows[tag] += edge_flow
873           
874               
875        return boundary_flows, total_boundary_inflow, total_boundary_outflow
876       
877
878    def compute_forcing_flows(self):
879        """Compute flows in and out of domain due to forcing terms.
880                       
881        Quantities computed are:
882               
883       
884           Total inflow through forcing terms
885           Total outflow through forcing terms
886           Current total volume in domain       
887
888        """
889
890        #FIXME(Ole): We need to separate what part of explicit_update was
891        # due to the normal flux calculations and what is due to forcing terms.
892       
893        pass
894                       
895       
896    def compute_total_volume(self):
897        """Compute total volume (m^3) of water in entire domain
898        """
899       
900        area = self.mesh.get_areas()
901        volume = 0.0
902       
903        stage = self.get_quantity('stage').get_values(location='centroids')
904        elevation = self.get_quantity('elevation').get_values(location='centroids')       
905        depth = stage-elevation
906       
907        return num.sum(depth*area)
908       
909       
910    def volumetric_balance_statistics(self):               
911        """Create volumetric balance report suitable for printing or logging.
912        """
913       
914        (boundary_flows, total_boundary_inflow,
915         total_boundary_outflow) = self.compute_boundary_flows() 
916       
917        s = '---------------------------\n'       
918        s += 'Volumetric balance report:\n'
919        s += '--------------------------\n'
920        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
921        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow       
922        s += 'Net boundary flow by tags [m^3/s]\n'
923        for tag in boundary_flows:
924            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
925       
926        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 
927        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
928       
929        # The go through explicit forcing update and record the rate of change for stage and
930        # record into forcing_inflow and forcing_outflow. Finally compute integral
931        # of depth to obtain total volume of domain.
932       
933        # FIXME(Ole): This part is not yet done.               
934       
935        return s       
936           
937################################################################################
938# End of class Shallow Water Domain
939################################################################################
940
941#-----------------
942# Flux computation
943#-----------------
944
945## @brief Compute fluxes and timestep suitable for all volumes in domain.
946# @param domain The domain to calculate fluxes for.
947def compute_fluxes(domain):
948    """Compute fluxes and timestep suitable for all volumes in domain.
949
950    Compute total flux for each conserved quantity using "flux_function"
951
952    Fluxes across each edge are scaled by edgelengths and summed up
953    Resulting flux is then scaled by area and stored in
954    explicit_update for each of the three conserved quantities
955    stage, xmomentum and ymomentum
956
957    The maximal allowable speed computed by the flux_function for each volume
958    is converted to a timestep that must not be exceeded. The minimum of
959    those is computed as the next overall timestep.
960
961    Post conditions:
962      domain.explicit_update is reset to computed flux values
963      domain.timestep is set to the largest step satisfying all volumes.
964
965    This wrapper calls the underlying C version of compute fluxes
966    """
967
968    import sys
969    from shallow_water_ext import compute_fluxes_ext_central \
970                                  as compute_fluxes_ext
971
972    N = len(domain)    # number_of_triangles
973
974    # Shortcuts
975    Stage = domain.quantities['stage']
976    Xmom = domain.quantities['xmomentum']
977    Ymom = domain.quantities['ymomentum']
978    Bed = domain.quantities['elevation']
979
980    timestep = float(sys.maxint)
981
982    flux_timestep = compute_fluxes_ext(timestep,
983                                       domain.epsilon,
984                                       domain.H0,
985                                       domain.g,
986                                       domain.neighbours,
987                                       domain.neighbour_edges,
988                                       domain.normals,
989                                       domain.edgelengths,
990                                       domain.radii,
991                                       domain.areas,
992                                       domain.tri_full_flag,
993                                       Stage.edge_values,
994                                       Xmom.edge_values,
995                                       Ymom.edge_values,
996                                       Bed.edge_values,
997                                       Stage.boundary_values,
998                                       Xmom.boundary_values,
999                                       Ymom.boundary_values,
1000                                       Stage.explicit_update,
1001                                       Xmom.explicit_update,
1002                                       Ymom.explicit_update,
1003                                       domain.already_computed_flux,
1004                                       domain.max_speed,
1005                                       int(domain.optimise_dry_cells))
1006
1007    domain.flux_timestep = flux_timestep
1008
1009################################################################################
1010# Module functions for gradient limiting
1011################################################################################
1012
1013##
1014# @brief Wrapper for C version of extrapolate_second_order_sw.
1015# @param domain The domain to operate on.
1016# @note MH090605 The following method belongs to the shallow_water domain class
1017#       see comments in the corresponding method in shallow_water_ext.c
1018def extrapolate_second_order_sw(domain):
1019    """Wrapper calling C version of extrapolate_second_order_sw"""
1020
1021    import sys
1022    from shallow_water_ext import extrapolate_second_order_sw as extrapol2
1023
1024    N = len(domain) # number_of_triangles
1025
1026    # Shortcuts
1027    Stage = domain.quantities['stage']
1028    Xmom = domain.quantities['xmomentum']
1029    Ymom = domain.quantities['ymomentum']
1030    Elevation = domain.quantities['elevation']
1031
1032    extrapol2(domain,
1033              domain.surrogate_neighbours,
1034              domain.number_of_boundaries,
1035              domain.centroid_coordinates,
1036              Stage.centroid_values,
1037              Xmom.centroid_values,
1038              Ymom.centroid_values,
1039              Elevation.centroid_values,
1040              domain.vertex_coordinates,
1041              Stage.vertex_values,
1042              Xmom.vertex_values,
1043              Ymom.vertex_values,
1044              Elevation.vertex_values,
1045              int(domain.optimise_dry_cells))
1046
1047##
1048# @brief Distribution from centroids to vertices specific to the SWW eqn.
1049# @param domain The domain to operate on.
1050def distribute_using_vertex_limiter(domain):
1051    """Distribution from centroids to vertices specific to the SWW equation.
1052
1053    It will ensure that h (w-z) is always non-negative even in the
1054    presence of steep bed-slopes by taking a weighted average between shallow
1055    and deep cases.
1056
1057    In addition, all conserved quantities get distributed as per either a
1058    constant (order==1) or a piecewise linear function (order==2).
1059
1060    FIXME: more explanation about removal of artificial variability etc
1061
1062    Precondition:
1063      All quantities defined at centroids and bed elevation defined at
1064      vertices.
1065
1066    Postcondition
1067      Conserved quantities defined at vertices
1068    """
1069
1070    # Remove very thin layers of water
1071    protect_against_infinitesimal_and_negative_heights(domain)
1072
1073    # Extrapolate all conserved quantities
1074    if domain.optimised_gradient_limiter:
1075        # MH090605 if second order,
1076        # perform the extrapolation and limiting on
1077        # all of the conserved quantities
1078
1079        if (domain._order_ == 1):
1080            for name in domain.conserved_quantities:
1081                Q = domain.quantities[name]
1082                Q.extrapolate_first_order()
1083        elif domain._order_ == 2:
1084            domain.extrapolate_second_order_sw()
1085        else:
1086            raise 'Unknown order'
1087    else:
1088        # Old code:
1089        for name in domain.conserved_quantities:
1090            Q = domain.quantities[name]
1091
1092            if domain._order_ == 1:
1093                Q.extrapolate_first_order()
1094            elif domain._order_ == 2:
1095                Q.extrapolate_second_order_and_limit_by_vertex()
1096            else:
1097                raise 'Unknown order'
1098
1099    # Take bed elevation into account when water heights are small
1100    balance_deep_and_shallow(domain)
1101
1102    # Compute edge values by interpolation
1103    for name in domain.conserved_quantities:
1104        Q = domain.quantities[name]
1105        Q.interpolate_from_vertices_to_edges()
1106
1107##
1108# @brief Distribution from centroids to edges specific to the SWW eqn.
1109# @param domain The domain to operate on.
1110def distribute_using_edge_limiter(domain):
1111    """Distribution from centroids to edges specific to the SWW eqn.
1112
1113    It will ensure that h (w-z) is always non-negative even in the
1114    presence of steep bed-slopes by taking a weighted average between shallow
1115    and deep cases.
1116
1117    In addition, all conserved quantities get distributed as per either a
1118    constant (order==1) or a piecewise linear function (order==2).
1119
1120
1121    Precondition:
1122      All quantities defined at centroids and bed elevation defined at
1123      vertices.
1124
1125    Postcondition
1126      Conserved quantities defined at vertices
1127    """
1128
1129    # Remove very thin layers of water
1130    protect_against_infinitesimal_and_negative_heights(domain)
1131
1132    for name in domain.conserved_quantities:
1133        Q = domain.quantities[name]
1134        if domain._order_ == 1:
1135            Q.extrapolate_first_order()
1136        elif domain._order_ == 2:
1137            Q.extrapolate_second_order_and_limit_by_edge()
1138        else:
1139            raise 'Unknown order'
1140
1141    balance_deep_and_shallow(domain)
1142
1143    # Compute edge values by interpolation
1144    for name in domain.conserved_quantities:
1145        Q = domain.quantities[name]
1146        Q.interpolate_from_vertices_to_edges()
1147
1148##
1149# @brief  Protect against infinitesimal heights and associated high velocities.
1150# @param domain The domain to operate on.
1151def protect_against_infinitesimal_and_negative_heights(domain):
1152    """Protect against infinitesimal heights and associated high velocities"""
1153
1154    from shallow_water_ext import protect
1155
1156    # Shortcuts
1157    wc = domain.quantities['stage'].centroid_values
1158    zc = domain.quantities['elevation'].centroid_values
1159    xmomc = domain.quantities['xmomentum'].centroid_values
1160    ymomc = domain.quantities['ymomentum'].centroid_values
1161
1162    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
1163            domain.epsilon, wc, zc, xmomc, ymomc)
1164
1165##
1166# @brief Wrapper for C function balance_deep_and_shallow_c().
1167# @param domain The domain to operate on.
1168def balance_deep_and_shallow(domain):
1169    """Compute linear combination between stage as computed by
1170    gradient-limiters limiting using w, and stage computed by
1171    gradient-limiters limiting using h (h-limiter).
1172    The former takes precedence when heights are large compared to the
1173    bed slope while the latter takes precedence when heights are
1174    relatively small.  Anything in between is computed as a balanced
1175    linear combination in order to avoid numerical disturbances which
1176    would otherwise appear as a result of hard switching between
1177    modes.
1178
1179    Wrapper for C implementation
1180    """
1181
1182    from shallow_water_ext import balance_deep_and_shallow \
1183                                  as balance_deep_and_shallow_c
1184
1185    # Shortcuts
1186    wc = domain.quantities['stage'].centroid_values
1187    zc = domain.quantities['elevation'].centroid_values
1188    wv = domain.quantities['stage'].vertex_values
1189    zv = domain.quantities['elevation'].vertex_values
1190
1191    # Momentums at centroids
1192    xmomc = domain.quantities['xmomentum'].centroid_values
1193    ymomc = domain.quantities['ymomentum'].centroid_values
1194
1195    # Momentums at vertices
1196    xmomv = domain.quantities['xmomentum'].vertex_values
1197    ymomv = domain.quantities['ymomentum'].vertex_values
1198
1199    balance_deep_and_shallow_c(domain,
1200                               wc, zc, wv, zv, wc,
1201                               xmomc, ymomc, xmomv, ymomv)
1202
1203
1204
1205################################################################################
1206# Standard forcing terms
1207################################################################################
1208
1209##
1210# @brief Apply gravitational pull in the presence of bed slope.
1211# @param domain The domain to apply gravity to.
1212# @note Wrapper for C function gravity_c().
1213def gravity(domain):
1214    """Apply gravitational pull in the presence of bed slope
1215    Wrapper calls underlying C implementation
1216    """
1217
1218    from shallow_water_ext import gravity as gravity_c
1219
1220    xmom_update = domain.quantities['xmomentum'].explicit_update
1221    ymom_update = domain.quantities['ymomentum'].explicit_update
1222
1223    stage = domain.quantities['stage']
1224    elevation = domain.quantities['elevation']
1225
1226    h = stage.centroid_values - elevation.centroid_values
1227    z = elevation.vertex_values
1228
1229    x = domain.get_vertex_coordinates()
1230    g = domain.g
1231
1232    gravity_c(g, h, z, x, xmom_update, ymom_update)    #, 1.0e-6)
1233
1234##
1235# @brief Apply friction to a surface (implicit).
1236# @param domain The domain to apply Manning friction to.
1237# @note Wrapper for C function manning_friction_c().
1238def manning_friction_implicit(domain):
1239    """Apply (Manning) friction to water momentum
1240    Wrapper for c version
1241    """
1242
1243    from shallow_water_ext import manning_friction_old
1244    from shallow_water_ext import manning_friction_new
1245
1246    xmom = domain.quantities['xmomentum']
1247    ymom = domain.quantities['ymomentum']
1248
1249    x = domain.get_vertex_coordinates()
1250   
1251    w = domain.quantities['stage'].centroid_values
1252    z = domain.quantities['elevation'].vertex_values
1253
1254    uh = xmom.centroid_values
1255    vh = ymom.centroid_values
1256    eta = domain.quantities['friction'].centroid_values
1257
1258    xmom_update = xmom.semi_implicit_update
1259    ymom_update = ymom.semi_implicit_update
1260
1261    N = len(domain)
1262    eps = domain.minimum_allowed_height
1263    g = domain.g
1264
1265    if domain.use_new_mannings:
1266        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
1267    else:
1268        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
1269   
1270
1271##
1272# @brief Apply friction to a surface (explicit).
1273# @param domain The domain to apply Manning friction to.
1274# @note Wrapper for C function manning_friction_c().
1275def manning_friction_explicit(domain):
1276    """Apply (Manning) friction to water momentum
1277    Wrapper for c version
1278    """
1279
1280    from shallow_water_ext import manning_friction_old
1281    from shallow_water_ext import manning_friction_new
1282
1283    xmom = domain.quantities['xmomentum']
1284    ymom = domain.quantities['ymomentum']
1285
1286    x = domain.get_vertex_coordinates()
1287   
1288    w = domain.quantities['stage'].centroid_values
1289    z = domain.quantities['elevation'].vertex_values
1290
1291    uh = xmom.centroid_values
1292    vh = ymom.centroid_values
1293    eta = domain.quantities['friction'].centroid_values
1294
1295    xmom_update = xmom.explicit_update
1296    ymom_update = ymom.explicit_update
1297
1298    N = len(domain)
1299    eps = domain.minimum_allowed_height
1300    g = domain.g
1301
1302
1303    if domain.use_new_mannings:
1304        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
1305    else:
1306        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
1307
1308
1309
1310# FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?)
1311##
1312# @brief Apply linear friction to a surface.
1313# @param domain The domain to apply Manning friction to.
1314# @note Is this still used (30 Oct 2007)?
1315def linear_friction(domain):
1316    """Apply linear friction to water momentum
1317
1318    Assumes quantity: 'linear_friction' to be present
1319    """
1320
1321    from math import sqrt
1322
1323    w = domain.quantities['stage'].centroid_values
1324    z = domain.quantities['elevation'].centroid_values
1325    h = w-z
1326
1327    uh = domain.quantities['xmomentum'].centroid_values
1328    vh = domain.quantities['ymomentum'].centroid_values
1329    tau = domain.quantities['linear_friction'].centroid_values
1330
1331    xmom_update = domain.quantities['xmomentum'].semi_implicit_update
1332    ymom_update = domain.quantities['ymomentum'].semi_implicit_update
1333
1334    N = len(domain) # number_of_triangles
1335    eps = domain.minimum_allowed_height
1336    g = domain.g #Not necessary? Why was this added?
1337
1338    for k in range(N):
1339        if tau[k] >= eps:
1340            if h[k] >= eps:
1341                S = -tau[k]/h[k]
1342
1343                #Update momentum
1344                xmom_update[k] += S*uh[k]
1345                ymom_update[k] += S*vh[k]
1346
1347def depth_dependent_friction(domain, default_friction,
1348                             surface_roughness_data,
1349                             verbose=False):
1350    """Returns an array of friction values for each wet element adjusted for depth.
1351
1352    Inputs:
1353        domain - computational domain object
1354        default_friction - depth independent bottom friction
1355        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
1356        friction region.
1357
1358    Outputs:
1359        wet_friction - Array that can be used directly to update friction as follows:
1360                       domain.set_quantity('friction', wet_friction)
1361
1362       
1363       
1364    """
1365
1366    import numpy as num
1367   
1368    # Create a temp array to store updated depth dependent friction for wet elements
1369    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
1370    N=len(domain)
1371    wet_friction    = num.zeros(N, num.float)
1372    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
1373   
1374   
1375    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
1376    # Recompute depth as vector 
1377    d = depth.get_values(location='centroids')
1378 
1379    # rebuild the 'friction' values adjusted for depth at this instant
1380    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
1381       
1382        # Get roughness data for each element
1383        n0 = float(surface_roughness_data[i,0])
1384        d1 = float(surface_roughness_data[i,1])
1385        n1 = float(surface_roughness_data[i,2])
1386        d2 = float(surface_roughness_data[i,3])
1387        n2 = float(surface_roughness_data[i,4])
1388       
1389       
1390        # Recompute friction values from depth for this element
1391               
1392        if d[i]   <= d1: 
1393            depth_dependent_friction = n1
1394        elif d[i] >= d2:
1395            depth_dependent_friction = n2
1396        else:
1397            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
1398           
1399        # check sanity of result
1400        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
1401            log.critical('%s >>>> WARNING: computed depth_dependent friction '
1402                         'out of range, ddf%f, n1=%f, n2=%f'
1403                         % (model_data.basename,
1404                            depth_dependent_friction, n1, n2))
1405       
1406        # update depth dependent friction  for that wet element
1407        wet_friction[i] = depth_dependent_friction
1408       
1409    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
1410   
1411    if verbose :
1412        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
1413        n_min=min(nvals)
1414        n_max=max(nvals)
1415       
1416        log.critical('         ++++ calculate_depth_dependent_friction - '
1417                     'Updated friction - range  %7.3f to %7.3f'
1418                     % (n_min, n_max))
1419   
1420    return wet_friction
1421
1422
1423
1424################################################################################
1425# Initialise module
1426################################################################################
1427
1428from anuga.utilities import compile
1429if compile.can_use_C_extension('shallow_water_ext.c'):
1430    # Underlying C implementations can be accessed
1431    from shallow_water_ext import assign_windfield_values
1432else:
1433    msg = 'C implementations could not be accessed by %s.\n ' % __file__
1434    msg += 'Make sure compile_all.py has been run as described in '
1435    msg += 'the ANUGA installation guide.'
1436    raise Exception, msg
1437
1438# Optimisation with psyco
1439from anuga.config import use_psyco
1440if use_psyco:
1441    try:
1442        import psyco
1443    except:
1444        import os
1445        if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']:
1446            pass
1447            #Psyco isn't supported on 64 bit systems, but it doesn't matter
1448        else:
1449            msg = ('WARNING: psyco (speedup) could not be imported, '
1450                   'you may want to consider installing it')
1451            log.critical(msg)
1452    else:
1453        psyco.bind(Domain.distribute_to_vertices_and_edges)
1454        psyco.bind(Domain.compute_fluxes)
1455
1456
1457if __name__ == "__main__":
1458    pass
Note: See TracBrowser for help on using the repository browser.