Ignore:
Timestamp:
Aug 30, 2004, 4:24:12 PM (21 years ago)
Author:
ole
Message:

Added C implementation of compute_fluxes

Location:
inundation/ga/storm_surge/pyvolution
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/config.py

    r234 r240  
    6161#use_extensions = False   #Do not use C-extensions
    6262
    63 #use_psyco = True  #Use psyco optimisations
    64 use_psyco = False  #Do not use psyco optimisations
     63use_psyco = True  #Use psyco optimisations
     64#use_psyco = False  #Do not use psyco optimisations
    6565
    6666

  • inundation/ga/storm_surge/pyvolution/doit.py

    r223 r240  
    1 
    2 
    31import os
    42
    5 
    6 
    7 os.system('python compile_all.py')
    8 
     3os.system('python compile.py')
    94os.system('python test_all.py')
    10 
    115os.system('python run_profile.py')
    12 
     6os.system('python show_balanced_limiters.py')

  • inundation/ga/storm_surge/pyvolution/domain.py

    r234 r240  
    7070
    7171        #Defaults
    72         from config import max_smallsteps, beta, newstyle
     72        from config import max_smallsteps, beta, newstyle, epsilon
    7373        self.beta = beta
     74        self.epsilon = epsilon       
    7475        self.newstyle = newstyle
    7576        self.default_order = 1
    … …  
    514515
    515516
    516 
    517 def compute_fluxes_c(domain, max_timestep):
    518     """Compute all fluxes and the timestep suitable for all volumes
    519     This version works directly with consecutive data structure
    520     and calls C-extension.
    521     Note that flux function is hardwired into C-extension.
    522     """
    523    
    524     import sys
    525 
    526     neighbours = Volume.neighbours
    527     neighbour_faces = Volume.neighbour_faces
    528     normals = Volume.normals
    529     flux = Volume.explicit_update   
    530 
    531     area = Volume.geometric[:,0]           
    532     radius = Volume.geometric[:,1]           
    533     edgelengths = Volume.geometric[:,2:]
    534 
    535 
    536     flux[:] = 0.0 #Reset stored fluxes to zero       
    537     from domain_ext import compute_fluxes
    538     timestep = compute_fluxes(domain.flux_function,
    539                               Volume.conserved_quantities_face0,
    540                               Volume.conserved_quantities_face1,
    541                               Volume.conserved_quantities_face2,
    542                               Volume.field_values_face0,
    543                               Volume.field_values_face1,
    544                               Volume.field_values_face2,
    545                               Boundary_value.conserved_quantities,
    546                               Boundary_value.field_values,
    547                               Vector.coordinates,
    548                               neighbours,
    549                               neighbour_faces,
    550                               normals,
    551                               flux, area, radius, edgelengths,
    552                               max_timestep)
    553 
    554     domain.timestep = timestep
    555 
    556 
    557 
    558    
    559517               
    560518##############################################
    … …  
    571529    #from python_versions import distribute_to_vertices_and_edges
    572530    #from python_versions import update_conserved_quantities
     531
     532
     533#Optimisation with psyco
     534from config import use_psyco
     535if use_psyco:
     536    try:
     537        import psyco
     538    except:
     539        msg = 'WARNING: psyco (speedup) could not import'+\
     540              ', you may want to consider installing it'
     541        print msg
     542    else:
     543        psyco.bind(Domain.distribute_to_vertices_and_edges)
     544        psyco.bind(Domain.update_boundary)
     545        psyco.bind(Domain.update_timestep)     #Not worth it
     546        psyco.bind(Domain.update_conserved_quantities)
     547
     548if __name__ == "__main__":
     549    pass

  • inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py

    r232 r240  
    22"""
    33
    4 from util import rotate
    54
    65class Boundary:

  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r234 r240  
    3030                                conserved_quantities, other_quantities)
    3131
    32         from config import minimum_allowed_height
     32        from config import minimum_allowed_height, g
    3333        self.minimum_allowed_height = minimum_allowed_height
     34        self.g = g
    3435
    3536        self.forcing_terms.append(gravity)
    … …  
    7778        distribute_to_vertices_and_edges(self)
    7879       
     80
     81
     82#Rotation of momentum vector
     83def rotate(q, normal, direction = 1):
     84    """Rotate the momentum component q (q[1], q[2])
     85    from x,y coordinates to coordinates based on normal vector.
     86
     87    If direction is negative the rotation is inverted.
     88   
     89    Input vector is preserved
     90
     91    This function is specific to the shallow water wave equation
     92    """
     93
     94    #FIXME: Needs to be tested
     95
     96    from Numeric import zeros, Float
     97   
     98    assert len(q) == 3,\
     99           'Vector of conserved quantities must have length 3'\
     100           'for 2D shallow water equation'
     101
     102    try:
     103        l = len(normal)
     104    except:
     105        raise 'Normal vector must be an Numeric array'
     106
     107    #FIXME: Put this test into C-extension as well
     108    assert l == 2, 'Normal vector must have 2 components'
     109
     110 
     111    n1 = normal[0]
     112    n2 = normal[1]   
     113   
     114    r = zeros(len(q), Float) #Rotated quantities
     115    r[0] = q[0]              #First quantity, height, is not rotated
     116
     117    if direction == -1:
     118        n2 = -n2
     119
     120
     121    r[1] =  n1*q[1] + n2*q[2]
     122    r[2] = -n2*q[1] + n1*q[2]
     123   
     124    return r
     125
     126
     127
    79128####################################
    80129# Flux computation       
    … …  
    93142
    94143    Bed elevations zl and zr.
    95     #FIXME: Remove those and pass in height directly
    96     #(unless we'll include a bed elevation discontinuity later)
    97    
    98144    """
    99145
    … …  
    101147    from math import sqrt
    102148    from Numeric import array
    103     from util import rotate
    104149
    105150    #Align momentums with x-axis
    … …  
    108153
    109154    z = (zl+zr)/2 #Take average of field values
    110     #FIXME: Maybe allow discontinuity later
    111155
    112156    w_left  = q_left[0]   #w=h+z
    … …  
    174218    Fluxes across each edge are scaled by edgelengths and summed up
    175219    Resulting flux is then scaled by area and stored in
    176     domain.explicit_update
     220    explicit_update for each of the three conserved quantities
     221    level, xmomentum and ymomentum   
    177222
    178223    The maximal allowable speed computed by the flux_function for each volume
    … …  
    187232    import sys
    188233    from Numeric import zeros, Float
    189     from config import max_timestep   
    190234
    191235    N = domain.number_of_elements
    192236   
    193     neighbours = domain.neighbours
    194     neighbour_edges = domain.neighbour_edges
    195     normals = domain.normals
    196 
    197     areas = domain.areas
    198     radii = domain.radii
    199     edgelengths = domain.edgelengths
    200        
    201     timestep = max_timestep #FIXME: Get rid of this
    202 
    203237    #Shortcuts
    204238    Level = domain.quantities['level']
    … …  
    220254
    221255    #Loop
     256    timestep = float(sys.maxint)   
    222257    for k in range(N):
    223         optimal_timestep = float(sys.maxint)
    224258
    225259        flux[:] = 0.  #Reset work array
    … …  
    230264
    231265            #Quantities at neighbour on nearest face
    232             n = neighbours[k,i]
     266            n = domain.neighbours[k,i]
    233267            if n < 0:
    234                 m = -n-1 #Convert neg flag to index
     268                m = -n-1 #Convert negative flag to index
    235269                qr = [level_bdry[m], xmom_bdry[m], ymom_bdry[m]]
    236270                zr = zl #Extend bed elevation to boundary
    237271            else:   
    238                 m = neighbour_edges[k,i]
     272                m = domain.neighbour_edges[k,i]
    239273                qr = [level[n, m], xmom[n, m], ymom[n, m]]
    240274                zr = bed[n, m]               
    … …  
    242276               
    243277            #Outward pointing normal vector   
    244             normal = normals[k, 2*i:2*i+2]
     278            normal = domain.normals[k, 2*i:2*i+2]
    245279
    246280            #Flux computation using provided function
    247281            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
    248             flux -= edgeflux * edgelengths[k,i]
    249            
     282            flux -= edgeflux * domain.edgelengths[k,i]
     283
    250284            #Update optimal_timestep
    251285            try:
    252                 optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
     286                timestep = min(timestep, domain.radii[k]/max_speed)
    253287            except ZeroDivisionError:
    254288                pass
    … …  
    256290        #Normalise by area and store for when all conserved
    257291        #quantities get updated
    258         flux /= areas[k]
     292        flux /= domain.areas[k]
    259293        Level.explicit_update[k] = flux[0]
    260294        Xmom.explicit_update[k] = flux[1]
    261295        Ymom.explicit_update[k] = flux[2]
    262296
    263         timestep = min(timestep, optimal_timestep)
    264 
    265297    domain.timestep = timestep   
    266298
    267        
     299
     300def compute_fluxes_c(domain):
     301    """Wrapper calling c version of compute fluxes
     302    """
     303
     304    import sys
     305    from Numeric import zeros, Float
     306
     307    N = domain.number_of_elements
     308   
     309    #Shortcuts
     310    Level = domain.quantities['level']
     311    Xmom = domain.quantities['xmomentum']
     312    Ymom = domain.quantities['ymomentum']
     313    Bed = domain.quantities['elevation']       
     314
     315    timestep = float(sys.maxint)   
     316    from shallow_water_ext import compute_fluxes
     317    domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
     318                                     domain.neighbours,
     319                                     domain.neighbour_edges,
     320                                     domain.normals,
     321                                     domain.edgelengths,                       
     322                                     domain.radii,
     323                                     domain.areas,
     324                                     Level.edge_values,
     325                                     Xmom.edge_values,
     326                                     Ymom.edge_values, 
     327                                     Bed.edge_values,   
     328                                     Level.boundary_values,
     329                                     Xmom.boundary_values,
     330                                     Ymom.boundary_values,
     331                                     Level.explicit_update,
     332                                     Xmom.explicit_update,
     333                                     Ymom.explicit_update)
     334       
    268335
    269336####################################
    … …  
    838905#Initialise module
    839906
     907
     908import compile
     909if compile.can_use_C_extension('shallow_water_ext.c'):
     910    #Replace python version with c implementations
     911    from shallow_water_ext import rotate
     912    compute_fluxes = compute_fluxes_c
     913    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
     914    #update_conserved_quantities = update_conserved_quantities_c   
     915
     916
     917#Optimisation with psyco
     918from config import use_psyco
     919if use_psyco:
     920    try:
     921        import psyco
     922    except:
     923        msg = 'WARNING: psyco (speedup) could not import'+\
     924              ', you may want to consider installing it'
     925        print msg
     926    else:
     927        psyco.bind(distribute_to_vertices_and_edges)
     928        psyco.bind(compute_fluxes_c) 
     929        #psyco.bind(update_boundary_values)
     930        #psyco.bind(Domain.update_timestep)     #Not worth it
     931        #psyco.bind(update_conserved_quantities)
     932
     933if __name__ == "__main__":
     934    pass

  • inundation/ga/storm_surge/pyvolution/test_shallow_water.py

    r232 r240  
    1717        pass
    1818        #print "  Tearing down"
     19
     20    def test_rotate(self):
     21        normal = array([0.0,-1.0])       
     22
     23        q = array([1.0,2.0,3.0])
     24     
     25        r = rotate(q, normal, direction = 1)
     26        assert r[0] == 1
     27        assert r[1] == -3
     28        assert r[2] == 2
     29
     30        w = rotate(r, normal, direction = -1)       
     31        assert allclose(w, q)
     32
    1933
    2034
    … …  
    189203       
    190204
    191     def test_compute_fluxes(self):
     205    def test_compute_fluxes0(self):
    192206        #Do a full triangle and check that fluxes cancel out for
    193207        #the constant level case

  • inundation/ga/storm_surge/pyvolution/test_util.py

    r205 r240  
    2727        assert zy == 0.0
    2828
    29 
    30     def test_rotate(self):
    31         normal = array([0.0,-1.0])       
    32 
    33         q = array([1.0,2.0,3.0])
    34      
    35         r = rotate(q, normal, direction = 1)
    36         assert r[0] == 1
    37         assert r[1] == -3
    38         assert r[2] == 2
    39 
    40         w = rotate(r, normal, direction = -1)       
    41         assert allclose(w, q)
    4229
    4330

  • inundation/ga/storm_surge/pyvolution/util.py

    r198 r240  
    8282
    8383
    84 def rotate_python(q, normal, direction = 1):
    85     """Rotate the momentum component q (q[1], q[2])
    86     from x,y coordinates to coordinates based on normal vector.
    87 
    88     If direction is negative the rotation is inverted.
    89    
    90     Input vector is preserved
    91 
    92     This function is specific to the shallow water wave equation
    93     """
    94 
    95     #FIXME: Needs to be tested
    96 
    97     from Numeric import zeros, Float
    98    
    99     assert len(q) == 3,\
    100            'Vector of conserved quantities must have length 3'\
    101            'for 2D shallow water equation'
    102 
    103     try:
    104         l = len(normal)
    105     except:
    106         raise 'Normal vector must be an Numeric array'
    107 
    108     #FIXME: Put this test into C-extension as well
    109     assert l == 2, 'Normal vector must have 2 components'
    110 
    111  
    112     n1 = normal[0]
    113     n2 = normal[1]   
    114    
    115     r = zeros(len(q), Float) #Rotated quantities
    116     r[0] = q[0]              #First quantity, height, is not rotated
    117 
    118     if direction == -1:
    119         n2 = -n2
    120 
    121 
    122     r[1] =  n1*q[1] + n2*q[2]
    123     r[2] = -n2*q[1] + n1*q[2]
    124    
    125     return r
    126 
    127 
    12884
    12985
    … …  
    13389import compile
    13490if compile.can_use_C_extension('util_ext.c'):
    135     from util_ext import gradient, rotate
     91    from util_ext import gradient
    13692else:
    13793    gradient = gradient_python
    138     rotate = rotate_python
    139 
    140 
    14194
    14295

  • inundation/ga/storm_surge/pyvolution/util_ext.c

    r205 r240  
    3636
    3737
    38 // Computational function for rotation
    39 int _rotate(double *q, double *r, double *normal, int direction) {
    40   /*Rotate the momentum component q (q[1], q[2])
    41     from x,y coordinates to coordinates based on normal vector.
    42    
    43     To rotate in opposite direction, call rotate with direction=1
    44     Result is returned in r
    45    
    46     q and r are assumed to have three components each*/
    47 
    48 
    49   double n0, n1, q1, q2;
    50  
    51   //Determine normal and direction
    52   n0 = normal[0];
    53   if (direction == 1) {
    54     n1 = normal[1];   
    55   } else {
    56     n1 = -normal[1];
    57   }
    58 
    59   //Shorthands
    60   q1 = q[1];  //uh momentum
    61   q2 = q[2];  //vh momentum
    62 
    63   //Rotate
    64   r[0] = q[0];
    65   r[1] = n0*q1 + n1*q2;
    66   r[2] = n0*q2 - n1*q1;
    67 
    68   return 0;
    69 } 
    70 
    71 
    72 
    7338// Gateway to Python
    7439PyObject *gradient(PyObject *self, PyObject *args) {
    … …  
    9661
    9762
    98 // Gateway to Python
    99 PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
    100   //
    101   // r = rotate(q, normal, direction=0)
    102   //
    103   // Where q is assumed to be a Float numeric array of length 3 and
    104   // normal a Float numeric array of length 2.
    105 
    106   //FIXME: Input checks and unit tests
    107  
    108   PyObject *Q, *Normal;
    109   PyArrayObject *q, *r, *normal;
    110  
    111   static char *argnames[] = {"q", "normal", "direction", NULL};
    112   int dimensions[1];
    113   int N, direction=0;
    114 
    115   // Convert Python arguments to C 
    116   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
    117                                    &Q, &Normal, &direction)) 
    118     return NULL; 
    119 
    120   //Input checks (convert sequenced into numeric arrays)
    121  
    122   q = (PyArrayObject *)
    123     PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
    124   normal = (PyArrayObject *)
    125     PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
    126  
    127   //Allocate space for return vector r (don't DECREF)
    128   N = q -> dimensions[0];
    129   dimensions[0] = N;
    130   r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    131 
    132   //Rotate
    133   _rotate((double *) q -> data,
    134           (double *) r -> data,
    135           (double *) normal -> data,     
    136           direction);     
    137 
    138   //Build result and DECREF r - returning r directly causes memory leak
    139   //result = Py_BuildValue("O", r);       
    140   //Py_DECREF(r);
    141  
    142   //return result;
    143   return PyArray_Return(r);
    144 }   
    145 
    146 
    147 // Method table for python module
    148 //static struct PyMethodDef MethodTable[] = {
    149 //  {"gradient", gradient, METH_VARARGS}, 
    150 //  {NULL, NULL}
    151 //};
    152 
    153 
    154 
    15563// Method table for python module
    15664static struct PyMethodDef MethodTable[] = {
    … …  
    16068   */
    16169 
    162   {"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"}, 
     70  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"}, 
    16371  {"gradient", gradient, METH_VARARGS, "Print out"},   
    16472  {NULL, NULL, 0, NULL}   /* sentinel */
Note: See TracChangeset for help on using the changeset viewer.