Changeset 4666
- Timestamp:
- Aug 8, 2007, 2:29:01 PM (18 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 4 edited
-
fit_interpolate/test_search_functions.py (modified) (1 diff)
-
shallow_water/benchmark_sww2dem.py (modified) (5 diffs)
-
utilities/quad.py (modified) (9 diffs)
-
utilities/test_quad.py (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/fit_interpolate/test_search_functions.py
r4655 r4666 182 182 mesh = Mesh(points, vertices) 183 183 184 # Don't do this, want to control the m x and mins184 # Don't do this, want to control the max and mins 185 185 #root = build_quadtree(mesh, max_points_per_cell=4) 186 186 187 #Initialise 188 Cell.initialise(mesh) 189 190 root = Cell(-3, 9, -3, 9, 187 188 root = Cell(-3, 9, -3, 9, mesh, 191 189 max_points_per_cell = 4) 192 190 #Insert indices of all vertices -
anuga_core/source/anuga/shallow_water/benchmark_sww2dem.py
r4664 r4666 22 22 import profile , pstats 23 23 import tempfile 24 import gc 25 24 26 from Scientific.IO.NetCDF import NetCDFFile 25 27 from Numeric import array … … 27 29 from anuga.fit_interpolate.interpolate import Interpolate 28 30 from anuga.fit_interpolate.fit import Fit 29 from anuga.pmesh.mesh import Mesh30 31 from anuga.geospatial_data.geospatial_data import Geospatial_data 31 32 from anuga.shallow_water import Domain … … 56 57 57 58 58 class Benchmarksww2dem:59 59 60 def 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 87 def build_sww(vert_rows, vert_columns, save): 60 88 """ 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 65 91 """ 66 67 def __init__(self):68 pass69 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 time92 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 79 105 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) 81 115 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 137 127 138 128 #------------------------------------------------------------- 139 129 if __name__ == "__main__": 140 b = Benchmarksww2dem()141 130 delimiter = ',' 142 131 143 132 ofile = 'lbm_resultsII.csv' 144 133 run_profile = False #True 145 size_list = [[4,5],[50,40]] 134 size_list = [[4,5],[5,40]] 135 #size_list = [[5,4]] 146 136 147 137 quantitiy_list = [['elevation'], … … 165 155 #sys.gettotalrefcount() 166 156 #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] 168 160 ,size[1] 169 161 ,output_quantities=quantitiy … … 172 164 print "time",time 173 165 print "mem", mem 166 gc.collect() 167 print "mem", mem 168 print " gc.get_objects() ", gc.get_objects() 174 169 #sys.gettotalrefcount() 175 170 fd.write(str(size[0]*size[1]) + delimiter + -
anuga_core/source/anuga/utilities/quad.py
r4655 r4666 27 27 """ 28 28 29 def __init__(self, southern, northern, western, eastern, 29 def __init__(self, southern, northern, western, eastern, mesh, 30 30 name = 'cell', 31 31 max_points_per_cell = 4): … … 39 39 self.western = round(western,5) 40 40 self.eastern = round(eastern,5) 41 self.mesh = mesh 41 42 42 43 # The points in this cell … … 63 64 cn = self.northern 64 65 cw = self.western 65 ce = self.eastern 66 ce = self.eastern 67 mesh = self.mesh 66 68 67 69 # 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')) 72 74 73 75 … … 140 142 if len(args) == 2: 141 143 point_id = int(args[1]) 142 x, y = self. __class__.mesh.get_nodes()[point_id]144 x, y = self.mesh.get_nodes()[point_id] 143 145 144 146 #print point_id, x, y … … 181 183 self.store(point) 182 184 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] 185 187 print "(" + str(x) + "," + str(y) + ")" 186 188 raise 'point not in region: %s' %str(point) … … 223 225 # print "verts", verts 224 226 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) 226 228 for k, _ in triangle_list: 227 229 if not triangles.has_key(k): 228 230 # 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) 233 235 triangles[k]=(tri, (n0, n1, n2)) 234 236 self.triangles = triangles.items() … … 418 420 419 421 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) 425 429 426 430 def build_quadtree(mesh, max_points_per_cell = 4): 427 431 """Build quad tree for mesh. 428 432 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. 430 435 """ 431 436 432 437 from Numeric import minimum, maximum 433 438 434 #Initialise435 Cell.initialise(mesh)436 439 437 440 #Make root cell … … 445 448 446 449 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 451 455 ymax += (ymax-ymin)/10 452 456 xmax += (xmax-xmin)/10 … … 462 466 463 467 #FIXME: Use mesh.filename if it exists 464 root = Cell(ymin, ymax, xmin, xmax, 468 root = Cell(ymin, ymax, xmin, xmax,mesh, 465 469 #name = .... 466 470 max_points_per_cell = max_points_per_cell) -
anuga_core/source/anuga/utilities/test_quad.py
r4655 r4666 12 12 13 13 def setUp(self): 14 self.cell = Cell(100, 140, 0, 40, 'cell')15 14 16 15 a = [3, 107] … … 30 29 mesh = Mesh(points, vertices) 31 30 self.mesh = mesh 32 Cell.initialise(mesh)31 self.cell = Cell(100, 140, 0, 40, mesh, 'cell') 33 32 34 33 def tearDown(self):
Note: See TracChangeset
for help on using the changeset viewer.
