Changeset 1855
- Timestamp:
- Sep 28, 2005, 6:00:46 PM (19 years ago)
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/pypar/mandelbrot_example/compile.py
r1852 r1855 111 111 include = os.path.join(os.path.join(sys.exec_prefix, 'include'), 112 112 'python'+sys.version[:3]) 113 headerfile = os.path.join(include, 'Python.h')114 113 114 115 # Check existence of Python.h 116 # 117 headerfile = include + '/Python.h' 115 118 try: 116 # unix style include directory117 119 open(headerfile, 'r') 118 120 except: 119 # windows style include directory 120 include = os.path.join(sys.exec_prefix, 'include') 121 headerfile = os.path.join(include, 'Python.h') 122 try: 123 open(headerfile, 'r') 124 except: 125 raise """Did not find Python header file. 126 Make sure files for Python C-extensions are installed.""" 121 raise """Did not find Python header file %s. 122 Make sure files for Python C-extensions are installed. 123 In debian linux, for example, you need to install a 124 package called something like python2.1-dev""" %headerfile 125 127 126 128 127 # Check filename(s) … … 147 146 # Compile 148 147 # 148 s = "%s -c %s -I%s -o %s.o" %(compiler, FN, include, root) 149 if os.name == 'posix' and os.uname()[4] == 'x86_64': 150 #Extra flags for 64 bit architectures 151 #s += ' -fPIC -m64' #gcc 152 s += ' -fPIC -tp amd64' #pgcc 153 149 154 150 s = "%s -c %s -I%s -o %s.o -Wall -O" %(compiler, FN, include, root)151 155 if verbose: 152 156 print s … … 161 165 162 166 # Make shared library (*.so) 167 s = "%s -%s %s -o %s.so" %(loader, sharedflag, object_files, root1) 168 169 if os.name == 'posix' and os.uname()[4] == 'x86_64': 170 #Extra flags for 64 bit architectures using Portland compilers 171 s += ' -mcmodel=medium' 163 172 164 s = "%s -%s %s -o %s.so -lm" %(loader, sharedflag, object_files, root1)165 173 if verbose: 166 174 print s -
inundation/pypar/mandelbrot_example/mandel_ext.c
r1852 r1855 1 / /Python - C extension for fast computation of the Mandelbrot iteration1 /* Python - C extension for fast computation of the Mandelbrot iteration 2 2 // 3 3 // To compile: … … 12 12 // k = calculate_point(c, kmax) 13 13 // 14 // Ole Nielsen, SUT 2003 14 // Ole Nielsen, SUT 2003 */ 15 15 16 16 #include "Python.h" 17 17 18 / / Computational function18 /* Computational function */ 19 19 int _calculate_point(Py_complex c, int kmax) { 20 20 int count; … … 37 37 38 38 39 / / Interface to Python39 /* Interface to Python */ 40 40 PyObject *calculate_point(PyObject *self, PyObject *args) { 41 PyComplexObject *C; / / Python Complex object42 Py_complex c; / / C Complex structure41 PyComplexObject *C; /* Python Complex object */ 42 Py_complex c; /* C Complex structure */ 43 43 int kmax, count; 44 44 45 / / Convert Python arguments to C45 /* Convert Python arguments to C */ 46 46 if (!PyArg_ParseTuple(args, "Oi", &C, &kmax)) 47 47 return NULL; 48 48 c = PyComplex_AsCComplex((PyObject*) C); 49 49 50 / / Call underlying routine50 /* Call underlying routine */ 51 51 count = _calculate_point(c, kmax); 52 52 53 / / Return results as a Python object53 /* Return results as a Python object */ 54 54 return Py_BuildValue("i", count); 55 55 } 56 56 57 57 58 / / Method table for python module58 /* Method table for python module */ 59 59 static struct PyMethodDef MethodTable[] = { 60 60 {"calculate_point", calculate_point, METH_VARARGS}, … … 63 63 64 64 65 / / Module initialisation65 /* Module initialisation */ 66 66 void initmandel_ext(void){ 67 67 Py_InitModule("mandel_ext", MethodTable); -
production/karratha_2005/create_mesh.py
r1837 r1855 11 11 12 12 13 14 def create_mesh(polygon, tags, resolution, filename = None): 15 """Create mesh from bounding polygon, tags for all segments and resolution. 16 17 Polygon is a list of points in latitude and longitudes - decimal degrees 18 19 20 Tags is a list of symbolic tags 21 22 Resolution is the maximal area 23 24 25 FIXME - FUTURE: 26 Use class Point_set 27 Take multiple polygons 28 29 Accept deg, min, sec, e.g. 30 [ [(-20,30,55), (116, 20, 00)], ...] 31 32 33 """ 34 35 from pmesh.mesh import * 36 from pyvolution.coordinate_transforms.redfearn import redfearn 37 38 39 #Convert to UTM 40 import project 41 refzone = project.refzone #FIXME 42 43 13 def convert_points_from_latlon_to_utm(polygon, refzone): 44 14 points = [] 45 15 for point in polygon: … … 48 18 points.append([easting, northing]) 49 19 50 polygon =points20 return points 51 21 52 53 mesh_origin = project.mesh_origin54 22 55 geo = Geo_reference(xllcorner = mesh_origin[1], #From dem 56 yllcorner = mesh_origin[2], 57 zone = refzone) 58 23 24 25 26 def create_region(polygon, tags, refzone): 27 """Create dictionary representation of region bounded by given polygon for use with pmesh 28 """ 29 59 30 60 31 # … … 77 48 #Create tags 78 49 #E.g. ['wall', 'wall', 'ocean', 'wall'] 79 segment_tags = [0]*N 80 for key in tags: 81 indices = tags[key] 82 for i in indices: 83 segment_tags[i] = key 50 51 if tags is not None: 52 segment_tags = [0]*N 53 for key in tags: 54 indices = tags[key] 55 for i in indices: 56 segment_tags[i] = key 84 57 85 dict['segment_tags'] = segment_tags58 dict['segment_tags'] = segment_tags 86 59 60 61 62 return dict 63 64 65 66 def create_mesh(bounding_polygon, boundary_tags, resolution, 67 filename = None, interior_regions = None): 68 """Create mesh from bounding polygon, tags for all segments and resolution. 69 70 Polygon is a list of points in latitude and longitudes - decimal degrees 71 72 73 Tags is a list of symbolic tags 74 75 Resolution is the maximal area 76 77 Interior_regions is a list of tuples consisting of (polygon, resolution) 78 79 This function does not allow segments to share points - use underlying 80 pmesh functionality for that 81 82 FIXME - FUTURE: 83 Use class Point_set 84 Take multiple polygons 85 86 Accept deg, min, sec, e.g. 87 [ [(-20,30,55), (116, 20, 00)], ...] 88 89 90 FIXME: This whole thing needs to be completely rethought 91 92 """ 93 94 from pmesh.mesh import * 95 from pyvolution.coordinate_transforms.redfearn import redfearn 96 from pyvolution.util import populate_polygon 97 98 99 #Make georef 100 #FIXME: Pass in geo or create automatically somehow 101 import project 102 refzone = project.refzone #FIXME 103 mesh_origin = project.mesh_origin 104 geo = Geo_reference(xllcorner = mesh_origin[1], #From dem 105 yllcorner = mesh_origin[2], 106 zone = refzone) 87 107 88 108 … … 92 112 93 113 m = Mesh(geo_reference=geo) 94 95 m.addVertsSegs(dict)96 114 97 115 98 #FIXME: Find centroid automatically 99 inner = m.addRegionEN(points[0][0]+1, points[0][1]+1) 116 117 #Convert to UTM 118 bounding_polygon = convert_points_from_latlon_to_utm(bounding_polygon, refzone) 119 120 121 #Do bounding polygon 122 region_dict = create_region(bounding_polygon, boundary_tags, refzone) 123 m.addVertsSegs(region_dict) 124 125 126 #Find one point inside region automatically 127 if interior_regions is not None: 128 excluded_polygons = [] 129 for P, res in interior_regions: 130 polygon = convert_points_from_latlon_to_utm(P, refzone) 131 excluded_polygons.append( polygon ) 132 else: 133 excluded_polygons = None 134 135 from Numeric import array 136 137 [inner_point] = populate_polygon(bounding_polygon, 1, exclude = excluded_polygons) 138 print 'I:', inner_point, resolution 139 inner = m.addRegionEN(inner_point[0], inner_point[1]) 100 140 inner.setMaxArea(resolution) 141 142 143 144 #Do interior regions 145 if interior_regions is not None: 146 for P, res in interior_regions: 147 polygon = convert_points_from_latlon_to_utm(P, refzone) 148 region_dict = create_region(polygon, None, refzone) 149 150 m.addVertsSegs(region_dict) 151 152 [inner_point] = populate_polygon(polygon, 1) 153 print 'I', inner_point, res 154 X = m.addRegionEN(inner_point[0], inner_point[1]) 155 X.setMaxArea(res) 156 157 101 158 102 159 m.generateMesh('pzq28.0za1000000a') -
production/karratha_2005/project.py
r1837 r1855 42 42 #east = degminsec2decimal_degrees(117,0,0) 43 43 44 south = degminsec2decimal_degrees(-20,5 0,0)44 south = degminsec2decimal_degrees(-20,55,0) 45 45 north = degminsec2decimal_degrees(-20,15,0) 46 west = degminsec2decimal_degrees(116, 20,0)46 west = degminsec2decimal_degrees(116,17,0) 47 47 east = degminsec2decimal_degrees(117,10,0) 48 48 … … 58 58 59 59 polygon = [p0, p1, p2, p3, p4, p5, p6, p7, p8] 60 refzone = 50 60 61 61 refzone = 50 62 63 #Interior regions 64 karratha_south = degminsec2decimal_degrees(-20,45,0) 65 karratha_north = degminsec2decimal_degrees(-20,37,0) 66 karratha_west = degminsec2decimal_degrees(116,45,0) 67 karratha_east = degminsec2decimal_degrees(116,55,0) 68 69 k0 = [karratha_south, karratha_west] 70 k1 = [karratha_south, karratha_east] 71 k2 = [karratha_north, karratha_east] 72 k3 = [karratha_north, karratha_west] 73 74 karratha_polygon = [k0, k1, k2, k3] 75 -
production/karratha_2005/run_karratha.py
r1845 r1855 10 10 Ole Nielsen, GA - 2005 11 11 """ 12 13 tide = 0.75 #HMWS estimate by Colin French, GA 14 12 15 13 16 import os … … 48 51 49 52 53 cache(ferret2sww, 54 (source_dir+project.boundary_basename, 55 project.boundary_basename), 56 {'verbose': True, 57 'minlat': south-1, 58 'maxlat': north+1, 59 'minlon': west-1, 60 'maxlon': east+1, 61 'origin': project.mesh_origin, 62 'mean_stage': tide, 63 'zscale': 1, 64 'fail_on_NaN': False, 65 'inverted_bathymetry': True}, 66 verbose = True) 67 #FIXME: Dependencies 50 68 51 #Create Triangular Mesh 52 from create_mesh import create_mesh 53 m = cache(create_mesh, 54 project.polygon, 55 {'tags': {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]}, 56 'resolution': 52000, 57 #'resolution': 150000, 58 'filename': project.meshname + '.msh'}, 59 verbose = True) 60 61 62 63 #tags = {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]}, 64 #resolution = 200000, 65 #filename = project.meshname + '.msh') 66 67 68 69 #ASSUME THAT SWW VERSION IS PRESENT 70 # 69 71 70 #ferret2sww(source_dir+project.boundary_basename, 72 71 # project.boundary_basename, … … 75 74 # minlon=west-1, maxlon=east+1, 76 75 # origin = project.mesh_origin, 76 # mean_stage = tide, 77 77 # zscale = 1, 78 78 # fail_on_NaN = False, 79 # inverted_bathymetry = True) ,79 # inverted_bathymetry = True) 80 80 81 82 #Create Triangular Mesh 83 from create_mesh import create_mesh 84 85 interior_regions = [[project.karratha_polygon, 10000]] 86 m = cache(create_mesh, 87 project.polygon, 88 {'boundary_tags': {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]}, 89 'resolution': 80000, 90 'filename': project.meshname + '.msh', 91 'interior_regions': interior_regions}, 92 verbose = True) 93 #verbose = True, 94 #evaluate = True ) 81 95 82 96 … … 92 106 domain.store = True 93 107 108 #domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] 109 domain.quantities_to_be_stored = ['stage'] 110 94 111 print 'Number of triangles = ', len(domain) 95 112 print 'The extent is ', domain.get_extent() … … 98 115 99 116 #Setup Initial Conditions 100 tide = 0101 117 domain.set_quantity('friction', 0) 102 118 domain.set_quantity('stage', tide) … … 112 128 113 129 Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True) 114 domain.starttime = 3000 #Obtained from MOST 130 #domain.starttime = 3000 #Obtained from MOST 131 domain.starttime = 0 #Obtained from MOST 115 132 116 133 Br = Reflective_boundary(domain) … … 126 143 127 144 #Run 145 #for t in domain.evolve(yieldstep = 600, finaltime = 15000): 146 # domain.write_time() 147 # domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage') 148 # 149 #for t in domain.evolve(yieldstep = 10, finaltime = 35000): 150 # domain.write_time() 151 # domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage') 152 128 153 for t in domain.evolve(yieldstep = 60, finaltime = 40000): 129 154 domain.write_time() 130 domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage') 155 domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')
Note: See TracChangeset
for help on using the changeset viewer.