Changeset 4666


Ignore:
Timestamp:
Aug 8, 2007, 2:29:01 PM (10 years ago)
Author:
duncan
Message:

unverified solution to ticket#176

Location:
anuga_core/source/anuga
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/test_search_functions.py

    r4655 r4666  
    182182        mesh = Mesh(points, vertices)
    183183
    184         # Don't do this, want to control the mx and mins
     184        # Don't do this, want to control the max and mins
    185185        #root = build_quadtree(mesh, max_points_per_cell=4)
    186186   
    187         #Initialise
    188         Cell.initialise(mesh)
    189 
    190         root = Cell(-3, 9, -3, 9,
     187
     188        root = Cell(-3, 9, -3, 9, mesh,
    191189                    max_points_per_cell = 4)
    192190        #Insert indices of all vertices
  • anuga_core/source/anuga/shallow_water/benchmark_sww2dem.py

    r4664 r4666  
    2222import profile , pstats
    2323import tempfile
     24import gc
     25
    2426from Scientific.IO.NetCDF import NetCDFFile
    2527from Numeric import array
     
    2729from anuga.fit_interpolate.interpolate import Interpolate
    2830from anuga.fit_interpolate.fit import Fit
    29 from anuga.pmesh.mesh import Mesh
    3031from anuga.geospatial_data.geospatial_data import Geospatial_data
    3132from anuga.shallow_water import Domain
     
    5657
    5758
    58 class Benchmarksww2dem:
    5959
     60def trial(vert_rows,
     61          vert_columns,
     62          output_quantities=['elevation'],
     63          save=False,
     64          verbose=False,
     65          run_profile=False):
     66    import time
     67 
     68    sww_file_name = build_sww(vert_rows, vert_columns, save=save)
     69       
     70    #Initial time and memory
     71    t0 = time.time()
     72    #m0 = None on windows
     73    m0 = mem_usage()
     74    #print "m0", m0
     75
     76    fileout = export_grid(sww_file_name, quantities=output_quantities)
     77       
     78    time_taken_sec = (time.time()-t0)
     79    m1 = mem_usage()
     80    #print "m1", m1
     81    if m0 is None or m1 is None:
     82        memory_used = None
     83    else:
     84        memory_used = (m1 - m0)
     85    return time_taken_sec, memory_used, m0, m1
     86
     87def build_sww(vert_rows, vert_columns, save):
    6088    """
    61 
    62     Note(DSG-DSG): If you are interested in benchmarking fitting, before
    63     and after blocking O:\1\dgray\before_blocking_subsandpit is before blocking
    64 
     89    Build an sww file we can read.
     90    verts_rows and vert_columns have to be greater that 1
    6591    """
    66 
    67     def __init__(self):
    68         pass
    69 
    70     def trial(self,
    71               vert_rows,
    72               vert_columns,
    73               output_quantities=['elevation'],
    74               save=False,
    75               verbose=False,
    76               run_profile=False):
    77  
    78         import time
     92    from anuga.pmesh.mesh import Mesh
     93    m = Mesh()
     94    m.build_grid(vert_rows, vert_columns)
     95    if save is True:
     96            m.export_mesh_file("aaaa.tsh")
     97    mesh_dict =  m.Mesh2IOTriangulationDict()
     98   
     99    sww_fileName = tempfile.mktemp(".sww" )
     100    # sww_fileName = "aa.sww"
     101    elevation = array(range(len(mesh_dict["vertices"])))
     102    stage = elevation
     103    ymomentum = elevation
     104    xmomentum = elevation
    79105       
    80         sww_file_name = self._build_sww(vert_rows, vert_columns, save=save)
     106    # NetCDF file definition
     107    fid = NetCDFFile(sww_fileName, 'w')
     108    sww = Write_sww()
     109    sww.header(fid, 0,
     110               len(mesh_dict['triangles']),
     111               len(mesh_dict["vertices"]))
     112    sww.triangulation(fid,
     113                      mesh_dict["vertices"], mesh_dict['triangles'],
     114                      elevation)
    81115       
    82         #Initial time and memory
    83         t0 = time.time()
    84         #m0 = None on windows
    85         m0 = mem_usage()
    86         #print "m0", m0
    87 
    88         fileout = export_grid(sww_file_name, quantities=output_quantities)
    89        
    90         time_taken_sec = (time.time()-t0)
    91         m1 = mem_usage()
    92         #print "m1", m1
    93         if m0 is None or m1 is None:
    94             memory_used = None
    95         else:
    96             memory_used = (m1 - m0)
    97         #print 'That took %.2f seconds' %time_taken_sec
    98         return time_taken_sec, memory_used, m0, m1
    99 
    100     def _build_sww(self, vert_rows, vert_columns, save):
    101         """
    102         Build an sww file we can read.
    103         verts_rows and vert_columns have to be greater that 1
    104         """
    105         m = Mesh()
    106         m.build_grid(vert_rows, vert_columns)
    107         if save is True:
    108             m.export_mesh_file("aaaa.tsh")
    109         mesh_dict =  m.Mesh2IOTriangulationDict()
    110        
    111         sww_fileName = tempfile.mktemp(".sww" )
    112         # sww_fileName = "aa.sww"
    113         elevation = array(range(len(mesh_dict["vertices"])))
    114         stage = elevation
    115         ymomentum = elevation
    116         xmomentum = elevation
    117        
    118         # NetCDF file definition
    119         fid = NetCDFFile(sww_fileName, 'w')
    120         sww = Write_sww()
    121         sww.header(fid, 0,
    122                    len(mesh_dict['triangles']),
    123                    len(mesh_dict["vertices"]))
    124         sww.triangulation(fid,
    125                    mesh_dict["vertices"], mesh_dict['triangles'],
    126                           elevation)
    127        
    128         for time_step in range(10):
    129             sww.quantities(fid,
    130                            time=time_step,
    131                            stage=stage,
    132                            xmomentum=xmomentum,
    133                            ymomentum=ymomentum)     
    134         #Close
    135         fid.close()
    136         return sww_fileName
     116    for time_step in range(10):
     117        sww.quantities(fid,
     118                       time=time_step,
     119                       stage=stage,
     120                       xmomentum=xmomentum,
     121                       ymomentum=ymomentum)     
     122    #Close
     123    fid.close()
     124    del mesh_dict
     125    del Mesh
     126    return sww_fileName
    137127   
    138128#-------------------------------------------------------------
    139129if __name__ == "__main__":
    140     b = Benchmarksww2dem()
    141130    delimiter = ','
    142131
    143132    ofile = 'lbm_resultsII.csv'
    144133    run_profile = False #True
    145     size_list = [[4,5],[50,40]]
     134    size_list = [[4,5],[5,40]]
     135    #size_list = [[5,4]]
    146136
    147137    quantitiy_list = [['elevation'],
     
    165155            #sys.gettotalrefcount()
    166156            #sys.getobjects()
    167             time, mem, m0, m1 = b.trial(size[0]
     157            gc.collect()
     158            #print " gc.get_objects() ", gc.get_objects()
     159            time, mem, m0, m1 = trial(size[0]
    168160                                        ,size[1]
    169161                                        ,output_quantities=quantitiy
     
    172164            print "time",time
    173165            print "mem", mem
     166            gc.collect()
     167            print "mem", mem
     168            print " gc.get_objects() ", gc.get_objects()
    174169            #sys.gettotalrefcount()
    175170            fd.write(str(size[0]*size[1]) + delimiter +
  • anuga_core/source/anuga/utilities/quad.py

    r4655 r4666  
    2727    """
    2828 
    29     def __init__(self, southern, northern, western, eastern,
     29    def __init__(self, southern, northern, western, eastern, mesh,
    3030                 name = 'cell',
    3131                 max_points_per_cell = 4):
     
    3939        self.western = round(western,5)   
    4040        self.eastern = round(eastern,5)
     41        self.mesh = mesh
    4142
    4243        # The points in this cell     
     
    6364        cn = self.northern
    6465        cw = self.western   
    65         ce = self.eastern
     66        ce = self.eastern   
     67        mesh = self.mesh
    6668
    6769        # create 4 child cells
    68         self.AddChild(Cell((cn+cs)/2,cn,cw,(cw+ce)/2,self.name+'_nw'))
    69         self.AddChild(Cell((cn+cs)/2,cn,(cw+ce)/2,ce,self.name+'_ne'))
    70         self.AddChild(Cell(cs,(cn+cs)/2,(cw+ce)/2,ce,self.name+'_se'))
    71         self.AddChild(Cell(cs,(cn+cs)/2,cw,(cw+ce)/2,self.name+'_sw'))
     70        self.AddChild(Cell((cn+cs)/2,cn,cw,(cw+ce)/2,mesh,self.name+'_nw'))
     71        self.AddChild(Cell((cn+cs)/2,cn,(cw+ce)/2,ce,mesh,self.name+'_ne'))
     72        self.AddChild(Cell(cs,(cn+cs)/2,(cw+ce)/2,ce,mesh,self.name+'_se'))
     73        self.AddChild(Cell(cs,(cn+cs)/2,cw,(cw+ce)/2,mesh,self.name+'_sw'))
    7274       
    7375 
     
    140142        if len(args) == 2:
    141143            point_id = int(args[1])
    142             x, y = self.__class__.mesh.get_nodes()[point_id]
     144            x, y = self.mesh.get_nodes()[point_id]
    143145
    144146            #print point_id, x, y
     
    181183                    self.store(point)
    182184                else:
    183                     x = self.__class__.mesh.coordinates[point][0]
    184                     y = self.__class__.mesh.coordinates[point][1]
     185                    x = self.mesh.coordinates[point][0]
     186                    y = self.mesh.coordinates[point][1]
    185187                    print "(" + str(x) + "," + str(y) + ")"
    186188                    raise 'point not in region: %s' %str(point)
     
    223225            # print "verts", verts
    224226            for vert in verts:
    225                 triangle_list = self.__class__.mesh.get_triangles_and_vertices_per_node(vert)
     227                triangle_list = self.mesh.get_triangles_and_vertices_per_node(vert)
    226228                for k, _ in triangle_list:
    227229                    if not triangles.has_key(k):
    228230                        # print 'k',k
    229                         tri = self.__class__.mesh.get_vertex_coordinates(k)
    230                         n0 = self.__class__.mesh.get_normal(k, 0)
    231                         n1 = self.__class__.mesh.get_normal(k, 1)
    232                         n2 = self.__class__.mesh.get_normal(k, 2)
     231                        tri = self.mesh.get_vertex_coordinates(k)
     232                        n0 = self.mesh.get_normal(k, 0)
     233                        n1 = self.mesh.get_normal(k, 1)
     234                        n2 = self.mesh.get_normal(k, 2)
    233235                        triangles[k]=(tri, (n0, n1, n2))
    234236            self.triangles = triangles.items()
     
    418420       
    419421
    420     #Class initialisation method       
    421     def initialise(cls, mesh):
    422         cls.mesh = mesh
    423 
    424     initialise = classmethod(initialise)
     422    #Class initialisation method
     423    # this is bad.  It adds a huge memory structure to the class.
     424    # When the instance is deleted the mesh hangs round (leaks).
     425    #def initialise(cls, mesh):
     426    #    cls.mesh = mesh
     427
     428    #initialise = classmethod(initialise)
    425429
    426430def build_quadtree(mesh, max_points_per_cell = 4):
    427431    """Build quad tree for mesh.
    428432
    429     All vertices in mesh are stored in quadtree and a reference to the root is returned.
     433    All vertices in mesh are stored in quadtree and a reference
     434    to the root is returned.
    430435    """
    431436
    432437    from Numeric import minimum, maximum
    433438
    434     #Initialise
    435     Cell.initialise(mesh)
    436439
    437440    #Make root cell
     
    445448
    446449   
    447     #Ensure boundary points are fully contained in region
    448     #It is a property of the cell structure that points on xmax or ymax of any given cell
    449     #belong to the neighbouring cell.
    450     #Hence, the root cell needs to be expanded slightly
     450    # Ensure boundary points are fully contained in region
     451    # It is a property of the cell structure that
     452    # points on xmax or ymax of any given cell
     453    # belong to the neighbouring cell.
     454    # Hence, the root cell needs to be expanded slightly
    451455    ymax += (ymax-ymin)/10
    452456    xmax += (xmax-xmin)/10
     
    462466   
    463467    #FIXME: Use mesh.filename if it exists
    464     root = Cell(ymin, ymax, xmin, xmax,
     468    root = Cell(ymin, ymax, xmin, xmax,mesh,
    465469                #name = ....
    466470                max_points_per_cell = max_points_per_cell)
  • anuga_core/source/anuga/utilities/test_quad.py

    r4655 r4666  
    1212
    1313    def setUp(self):
    14         self.cell = Cell(100, 140, 0, 40, 'cell')
    1514
    1615        a = [3, 107]
     
    3029        mesh = Mesh(points, vertices)
    3130        self.mesh = mesh
    32         Cell.initialise(mesh)
     31        self.cell = Cell(100, 140, 0, 40, mesh, 'cell')
    3332
    3433    def tearDown(self):
Note: See TracChangeset for help on using the changeset viewer.