Changeset 7452


Ignore:
Timestamp:
Sep 1, 2009, 10:29:19 AM (11 years ago)
Author:
steve
Message:

Committing latest parallel code

Location:
anuga_core/source
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/interface.py

    r7326 r7452  
    1818
    1919from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Time_boundary
     20from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Transmissive_boundary
     21
    2022from anuga.abstract_2d_finite_volumes.util import file_function
     23
     24from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    2125
    2226from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
     
    2630from anuga.utilities.polygon import Polygon_function
    2731
     32from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
     33
     34
     35#-----------------------------
     36# rectangular domains
     37#-----------------------------
     38def rectangular_cross_domain(*args, **kwargs):
     39    points, vertices, boundary = rectangular_cross(*args, **kwargs)
     40    return Domain(points, vertices, boundary)
     41   
     42#----------------------------
     43# Create domain from file
     44#----------------------------
     45def create_domain_from_file(file):
     46    return pmesh_to_domain_instance(file,Domain)
    2847
    2948#---------------------------
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r7352 r7452  
    557557    # @note 'polyline' may contain multiple sections allowing complex shapes.
    558558    # @note Assume absolute UTM coordinates.
    559     def new_get_energy_through_cross_section(self, polyline,
     559    def get_energy_through_cross_section(self, polyline,
    560560                                         kind='total',
    561561                                         verbose=False):
     
    588588        cross_section = Cross_section(self, polyline, verbose)
    589589
    590         return cross_section.get_energy_through_cross_section()
     590        return cross_section.get_energy_through_cross_section(kind)
    591591
    592592
     
    598598    # @note 'polyline' may contain multiple sections allowing complex shapes.
    599599    # @note Assume absolute UTM coordinates.
    600     def get_energy_through_cross_section(self, polyline,
     600    def old_get_energy_through_cross_section(self, polyline,
    601601                                         kind='total',
    602602                                         verbose=False):
     
    27862786        self.midpoints = ensure_geospatial(self.midpoints, self.domain.geo_reference)
    27872787
     2788    ##
     2789    # @brief set verbose mode
     2790    def set_verbose(self,verbose=True):
     2791        """Set verbose mode true or flase
     2792        """
     2793
     2794        self.verbose=verbose
    27882795
    27892796    ##
     
    28202827    ##
    28212828    # @brief calculate current energy flow through cross section
    2822     def get_energy_through_cross_section(self):
     2829    def get_energy_through_cross_section(self, kind='total'):
    28232830        """Obtain average energy head [m] across specified cross section.
    28242831
     
    28372844        """
    28382845
     2846        from anuga.config import g, epsilon, velocity_protection as h0
     2847       
    28392848        # Get interpolated values
    28402849        stage = self.domain.get_quantity('stage')
     
    28432852        ymomentum = self.domain.get_quantity('ymomentum')
    28442853
    2845         w = stage.get_values(interpolation_points=midpoints, use_cache=True)
    2846         z = elevation.get_values(interpolation_points=midpoints, use_cache=True)
    2847         uh = xmomentum.get_values(interpolation_points=midpoints,
     2854        w = stage.get_values(interpolation_points=self.midpoints, use_cache=True)
     2855        z = elevation.get_values(interpolation_points=self.midpoints, use_cache=True)
     2856        uh = xmomentum.get_values(interpolation_points=self.midpoints,
    28482857                                  use_cache=True)
    2849         vh = ymomentum.get_values(interpolation_points=midpoints,
     2858        vh = ymomentum.get_values(interpolation_points=self.midpoints,
    28502859                                  use_cache=True)
    28512860        h = w-z                # Depth
     
    28532862        # Compute total length of polyline for use with weighted averages
    28542863        total_line_length = 0.0
    2855         for segment in segments:
     2864        for segment in self.segments:
    28562865            total_line_length += segment.length
    28572866
     
    28792888
    28802889            # Add to weighted average
    2881             weigth = segments[i].length/total_line_length
     2890            weigth = self.segments[i].length/total_line_length
    28822891            average_energy += segment_energy*weigth
    28832892
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r7398 r7452  
    16441644                                                             verbose=False)
    16451645                assert num.allclose(Et, w + 0.5*u*u/g)
     1646
     1647
     1648    def test_cross_section_class(self):
     1649        """test_cross_section_class(self):
     1650
     1651        Test that the total and specific energy through a cross section can be
     1652        correctly obtained at run-time from the ANUGA cross section class.
     1653
     1654        This test creates a flat bed with a known flow through it, creates a cross
     1655        section and tests that the correct flow and energies are calculated
     1656
     1657        The specifics are
     1658        e = -1 m
     1659        u = 2 m/s
     1660        h = 2 m
     1661        w = 3 m (width of channel)
     1662
     1663        q = u*h*w = 12 m^3/s
     1664
     1665        This run tries it with georeferencing and with elevation = -1
     1666        """
     1667
     1668        import time, os
     1669        from Scientific.IO.NetCDF import NetCDFFile
     1670        from mesh_factory import rectangular
     1671
     1672        # Create basic mesh (20m x 3m)
     1673        width = 3
     1674        length = 20
     1675        t_end = 1
     1676        points, vertices, boundary = rectangular(length, width, length, width)
     1677
     1678        # Create shallow water domain
     1679        domain = Domain(points, vertices, boundary,
     1680                        geo_reference=Geo_reference(56, 308500, 6189000))
     1681
     1682        domain.default_order = 2
     1683        domain.set_quantities_to_be_stored(None)
     1684
     1685        e = -1.0
     1686        w = 1.0
     1687        h = w-e
     1688        u = 2.0
     1689        uh = u*h
     1690
     1691        Br = Reflective_boundary(domain)       # Side walls
     1692        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
     1693
     1694        # Initial conditions
     1695        domain.set_quantity('elevation', e)
     1696        domain.set_quantity('stage', w)
     1697        domain.set_quantity('xmomentum', uh)
     1698        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     1699
     1700        # Interpolation points down the middle
     1701        I = [[0, width/2.],
     1702             [length/2., width/2.],
     1703             [length, width/2.]]
     1704        interpolation_points = domain.geo_reference.get_absolute(I)
     1705
     1706        # Shortcuts to quantites
     1707        stage = domain.get_quantity('stage')
     1708        xmomentum = domain.get_quantity('xmomentum')
     1709        ymomentum = domain.get_quantity('ymomentum')
     1710
     1711
     1712        # Create some cross sections
     1713        cross_sections = []
     1714        for i in range(5):
     1715            x = length/2. + i*0.23674563    # Arbitrary
     1716            polyline = [[x, 0], [x, width]]
     1717
     1718            polyline = domain.geo_reference.get_absolute(polyline)
     1719       
     1720            cross_sections.append(Cross_section(domain,polyline))
     1721
     1722
     1723 
     1724        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
     1725            # Check that quantities are they should be in the interior
     1726            w_t = stage.get_values(interpolation_points)
     1727            uh_t = xmomentum.get_values(interpolation_points)
     1728            vh_t = ymomentum.get_values(interpolation_points)
     1729
     1730            assert num.allclose(w_t, w)
     1731            assert num.allclose(uh_t, uh)
     1732            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
     1733
     1734
     1735            # Check flows and energies through the middle
     1736            for cross_section in cross_sections:
     1737               
     1738                Q = cross_section.get_flow_through_cross_section()
     1739
     1740                assert num.allclose(Q, uh*width)
     1741
     1742                Es = cross_section.get_energy_through_cross_section(kind='specific')
     1743
     1744                assert num.allclose(Es, h + 0.5*u*u/g)
     1745
     1746                Et = cross_section.get_energy_through_cross_section(kind='total')
     1747               
     1748                assert num.allclose(Et, w + 0.5*u*u/g)
     1749
     1750
     1751
     1752
     1753
    16461754
    16471755    def test_another_runup_example(self):
  • anuga_core/source/anuga/visualiser/offline.py

    r4325 r7452  
    1 from Numeric import array, Float, ravel, zeros
     1#from Numeric import array, Float, ravel, zeros
     2import numpy as num
    23from Scientific.IO.NetCDF import NetCDFFile
    34from Tkinter import Button, E, Tk, W, Label, StringVar, Scale, HORIZONTAL
     
    8081        else:
    8182            N_vert = len(fin.variables[quantityName])
    82         x = ravel(array(fin.variables['x'], Float))
    83         y = ravel(array(fin.variables['y'], Float))
     83        x = num.ravel(num.array(fin.variables['x'], num.float))
     84        y = num.ravel(num.array(fin.variables['y'], num.float))
    8485        if dynamic is True:
    85             q = array(fin.variables[quantityName][frameNumber], Float)
    86         else:
    87             q = ravel(array(fin.variables[quantityName], Float))
     86            q = num.array(fin.variables[quantityName][frameNumber], num.float)
     87        else:
     88            q = num.ravel(num.array(fin.variables[quantityName], num.float))
    8889
    8990        q *= self.height_zScales[quantityName]
     
    123124        for q in filter(lambda n:n != 'x' and n != 'y' and n != 'z' and n != 'time' and n != 'volumes', fin.variables.keys()):
    124125            if len(fin.variables[q].shape) == 1: # Not a time-varying quantity
    125                 quantities[q] = ravel(array(fin.variables[q], Float))
     126                quantities[q] = num.ravel(num.array(fin.variables[q], num.float))
    126127            else: # Time-varying, get the current timestep data
    127                 quantities[q] = array(fin.variables[q][self.frameNumber], Float)
     128                quantities[q] = num.array(fin.variables[q][self.frameNumber], num.float)
    128129        fin.close()
    129130        return quantities
  • anuga_core/source/anuga/visualiser/realtime.py

    r6113 r7452  
    1 from Numeric import Float, zeros, shape
     1#from  import Float, zeros, shape
     2import numpy as num
    23from Tkinter import Button, E, Tk, W
    34from threading import Event
     
    4748
    4849        # Also build vert_index - a list of the x & y values of each vertex
    49         self.vert_index = zeros((N_vert,2), Float)
     50        self.vert_index = num.zeros((N_vert,2), num.float)
    5051        for n in range(N_tri):
    5152            self.vtk_cells.InsertNextCell(3)
     
    5657    def update_height_quantity(self, quantityName, dynamic=True):
    5758        N_vert = len(self.source.get_vertex_coordinates())
    58         qty_index = zeros(N_vert, Float)
     59        qty_index = num.zeros(N_vert, num.float)
    5960        triangles = self.source.get_triangles()
    6061        vertex_values, _ = self.source.get_quantity(quantityName).get_vertex_values(xy=False, smooth=False)
  • anuga_core/source/anuga_parallel/test_parallel_sw_flow.py

    r7449 r7452  
    1515import os
    1616import sys
    17 import pypar
     17#import pypar
    1818import numpy as num
    19 
    20 from anuga.utilities.numerical_tools import ensure_numeric
    21 from anuga.utilities.polygon import is_inside_polygon
    2219
    2320from anuga.interface import Domain
     
    3027
    3128
    32 from anuga_parallel.interface import distribute, myid, numprocs, send, receive
     29from anuga_parallel.interface import distribute, myid, numprocs, send, receive, barrier
    3330
    3431#--------------------------------------------------------------------------
     
    3835finaltime = 6.0
    3936nprocs = 4
    40 N = 29
    41 M = 29
    42 verbose = True
     37N = 25
     38M = 25
     39verbose =  True
    4340
    4441#---------------------------------
     
    5148# Setup Test
    5249##########################################################################
    53 def evolution_test(parallel=False, G = None):
     50def evolution_test(parallel=False, G = None, seq_interpolation_points=None):
    5451
    5552    #--------------------------------------------------------------------------
     
    6562    #--------------------------------------------------------------------------
    6663    if parallel:
    67         if myid == 0: print 'DISTRIBUTING PARALLEL DOMAIN'
     64        if myid == 0 and verbose : print 'DISTRIBUTING PARALLEL DOMAIN'
    6865        domain = distribute(domain, verbose=False)
    6966
     
    9188    #------------------------------------------------------------------------------
    9289    # Find which sub_domain in which the interpolation points are located
    93     #------------------------------------------------------------------------------
    94     interpolation_points = [[0.4,0.51], [0.6,0.51], [0.8,0.51], [0.9,0.51]]
     90    #
     91    # Sometimes the interpolation points sit exactly
     92    # between to centroids, so in the parallel run we
     93    # reset the interpolation points to the centroids
     94    # found in the sequential run
     95    #------------------------------------------------------------------------------
     96    interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
     97
     98    if parallel:
     99        interpolation_points = seq_interpolation_points
     100
     101
    95102    gauge_values = []
    96103    tri_ids = []
     
    108115            tri_ids.append(-2)
    109116
     117
    110118    if verbose: print 'P%d has points = %s' %(myid, tri_ids)
    111119
    112 
    113 
     120    if not parallel:
     121        c_coord = domain.get_centroid_coordinates()
     122        interpolation_points = []
     123        for id in tri_ids:
     124            if id<1:
     125                print 'ERROR: All interpolation points be within the sequential domain!'
     126            interpolation_points.append(c_coord[id,:])
     127           
    114128    #------------------------------------------------------------------------------
    115129    # Evolve system through time
     
    152166    assert_(success)
    153167   
    154     return G
     168    return G, interpolation_points
    155169
    156170# Test an nprocs-way run of the shallow water equations
     
    182196        # array G
    183197        #------------------------------------------
    184         pypar.barrier()
     198        barrier()
    185199        if myid == 0 and verbose: print 'SEQUENTIAL START'
    186200
    187         G = evolution_test(parallel=False)
     201        G , interpolation_points = evolution_test(parallel=False)
    188202        G = num.array(G,num.float)
    189203
    190         pypar.barrier()
     204        barrier()
    191205       
    192206        #------------------------------------------
     
    196210        if myid ==0 and verbose: print 'PARALLEL START'
    197211
    198         evolution_test(parallel=True, G=G)
     212        evolution_test(parallel=True, G=G, seq_interpolation_points = interpolation_points)
    199213       
    200214
  • anuga_core/source/pymetis/pymetis/metis.c

    r3414 r7452  
    11#include <Python.h>
    2 #include <Numeric/arrayobject.h>
     2#include <numpy/arrayobject.h>
    33
    44/* This must be the same as the metis idxtype */
     
    5050  int edgecut;
    5151  int numflag = 0; // The metis routine requires an int * for numflag.
    52   int dims[1]; // PyArray_FromDimsAndData needs an int[] of array sizes.
     52  npy_intp dims[1]; // PyArray_SimpleNewFromData needs an npy_intp[] of array sizes.
    5353
    5454  PyObject * elements;
     
    104104
    105105  dims[0] = ne;
    106   epart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)epart);
     106  epart_pyarr = (PyArrayObject *)PyArray_SimpleNewFromData(1, dims, PyArray_INT, (void *)epart);
     107  //epart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)epart);
    107108  dims[0] = nn;
    108   npart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)npart);
     109  npart_pyarr = (PyArrayObject *)PyArray_SimpleNewFromData(1, dims, PyArray_INT, (void *)npart);
     110  //npart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)npart);
    109111
    110112 
  • anuga_core/source/pymetis/test_metis.py

    r5846 r7452  
    33from sys import path
    44
    5 from Numeric import array, allclose
     5from numpy import array, allclose
    66
    77import unittest
Note: See TracChangeset for help on using the changeset viewer.