Changeset 1855


Ignore:
Timestamp:
Sep 28, 2005, 6:00:46 PM (19 years ago)
Author:
ole
Message:

Better interface to meshing for karratha + some pypar work

Files:
5 edited

Legend:

Unmodified
Added
Removed
  • inundation/pypar/mandelbrot_example/compile.py

    r1852 r1855  
    111111  include = os.path.join(os.path.join(sys.exec_prefix, 'include'),
    112112                     'python'+sys.version[:3])
    113   headerfile = os.path.join(include, 'Python.h')
    114113 
     114 
     115  # Check existence of Python.h
     116  #
     117  headerfile = include + '/Python.h'
    115118  try:
    116     # unix style include directory
    117119    open(headerfile, 'r')
    118120  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 
    127126 
    128127  # Check filename(s)
     
    147146    # Compile
    148147    #
     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     
    149154   
    150     s = "%s -c %s -I%s -o %s.o -Wall -O" %(compiler, FN, include, root)
    151155    if verbose:
    152156      print s
     
    161165 
    162166  # 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'
    163172 
    164   s = "%s -%s %s -o %s.so -lm" %(loader, sharedflag, object_files, root1)
    165173  if verbose:
    166174    print s
  • inundation/pypar/mandelbrot_example/mandel_ext.c

    r1852 r1855  
    1 // Python - C extension for fast computation of the Mandelbrot iteration
     1/* Python - C extension for fast computation of the Mandelbrot iteration
    22//
    33// To compile:
     
    1212// k = calculate_point(c, kmax)
    1313//
    14 // Ole Nielsen, SUT 2003
     14// Ole Nielsen, SUT 2003 */
    1515       
    1616#include "Python.h"
    1717
    18 // Computational function
     18/* Computational function */
    1919int _calculate_point(Py_complex c, int kmax) {
    2020  int count;
     
    3737
    3838
    39 // Interface to Python
     39/* Interface to Python */
    4040PyObject *calculate_point(PyObject *self, PyObject *args) {
    41   PyComplexObject *C;  // Python Complex object
    42   Py_complex c;        // C Complex structure
     41  PyComplexObject *C;  /* Python Complex object */
     42  Py_complex c;        /* C Complex structure */
    4343  int kmax, count;
    4444
    45   // Convert Python arguments to C 
     45  /* Convert Python arguments to C  */
    4646  if (!PyArg_ParseTuple(args, "Oi", &C, &kmax))
    4747    return NULL;
    4848  c = PyComplex_AsCComplex((PyObject*) C);   
    4949 
    50   // Call underlying routine
     50  /* Call underlying routine */
    5151  count = _calculate_point(c, kmax);
    5252
    53   // Return results as a Python object
     53  /* Return results as a Python object */
    5454  return Py_BuildValue("i", count);
    5555}
    5656
    5757
    58 // Method table for python module
     58/* Method table for python module */
    5959static struct PyMethodDef MethodTable[] = {
    6060  {"calculate_point", calculate_point, METH_VARARGS}, 
     
    6363
    6464
    65 // Module initialisation   
     65/* Module initialisation  */
    6666void initmandel_ext(void){
    6767  Py_InitModule("mandel_ext", MethodTable);
  • production/karratha_2005/create_mesh.py

    r1837 r1855  
    1111
    1212
    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 
     13def convert_points_from_latlon_to_utm(polygon, refzone):
    4414    points = []
    4515    for point in polygon:
     
    4818        points.append([easting, northing])
    4919
    50     polygon = points   
     20    return points   
    5121
    52    
    53     mesh_origin = project.mesh_origin
    5422
    55     geo = Geo_reference(xllcorner = mesh_origin[1],  #From dem
    56                         yllcorner = mesh_origin[2],
    57                         zone = refzone)
    58    
     23
     24
     25
     26def create_region(polygon, tags, refzone):
     27    """Create dictionary representation of region bounded by given polygon for use with pmesh
     28    """
     29
    5930
    6031    #
     
    7748    #Create tags
    7849    #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
    8457   
    85     dict['segment_tags'] = segment_tags
     58        dict['segment_tags'] = segment_tags
    8659
     60
     61
     62    return dict
     63
     64
     65
     66def 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)
    87107
    88108
     
    92112   
    93113    m = Mesh(geo_reference=geo)
    94        
    95     m.addVertsSegs(dict)
    96114
    97115
    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])
    100140    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
    101158   
    102159    m.generateMesh('pzq28.0za1000000a')
  • production/karratha_2005/project.py

    r1837 r1855  
    4242#east = degminsec2decimal_degrees(117,0,0)
    4343
    44 south = degminsec2decimal_degrees(-20,50,0)
     44south = degminsec2decimal_degrees(-20,55,0)
    4545north = degminsec2decimal_degrees(-20,15,0)
    46 west = degminsec2decimal_degrees(116,20,0)
     46west = degminsec2decimal_degrees(116,17,0)
    4747east = degminsec2decimal_degrees(117,10,0)
    4848
     
    5858   
    5959polygon = [p0, p1, p2, p3, p4, p5, p6, p7, p8]
     60refzone = 50
    6061
    61 refzone = 50
     62
     63#Interior regions
     64karratha_south = degminsec2decimal_degrees(-20,45,0)
     65karratha_north = degminsec2decimal_degrees(-20,37,0)
     66karratha_west = degminsec2decimal_degrees(116,45,0)
     67karratha_east = degminsec2decimal_degrees(116,55,0)
     68
     69k0 = [karratha_south, karratha_west]
     70k1 = [karratha_south, karratha_east]
     71k2 = [karratha_north, karratha_east]
     72k3 = [karratha_north, karratha_west]   
     73
     74karratha_polygon = [k0, k1, k2, k3]
     75
  • production/karratha_2005/run_karratha.py

    r1845 r1855  
    1010Ole Nielsen, GA - 2005
    1111"""
     12
     13tide = 0.75   #HMWS estimate by Colin French, GA
     14
    1215
    1316import os
     
    4851
    4952
     53cache(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
    5068
    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   
    7170#ferret2sww(source_dir+project.boundary_basename,
    7271#           project.boundary_basename,
     
    7574#           minlon=west-1, maxlon=east+1,
    7675#           origin = project.mesh_origin,
     76#           mean_stage = tide,
    7777#           zscale = 1,
    7878#           fail_on_NaN = False,
    79 #           inverted_bathymetry = True),
     79#           inverted_bathymetry = True)
    8080
     81
     82#Create Triangular Mesh
     83from create_mesh import create_mesh
     84
     85interior_regions = [[project.karratha_polygon, 10000]]
     86m = 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 )
    8195
    8296
     
    92106domain.store = True
    93107
     108#domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
     109domain.quantities_to_be_stored = ['stage']
     110       
    94111print 'Number of triangles = ', len(domain)
    95112print 'The extent is ', domain.get_extent()
     
    98115
    99116#Setup Initial Conditions
    100 tide = 0
    101117domain.set_quantity('friction', 0)
    102118domain.set_quantity('stage', tide)
     
    112128
    113129Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True)
    114 domain.starttime = 3000  #Obtained from MOST
     130#domain.starttime = 3000  #Obtained from MOST
     131domain.starttime = 0  #Obtained from MOST
    115132
    116133Br = Reflective_boundary(domain)
     
    126143
    127144#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
    128153for t in domain.evolve(yieldstep = 60, finaltime = 40000):
    129154    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.