Changeset 7823


Ignore:
Timestamp:
Jun 11, 2010, 6:14:08 PM (10 years ago)
Author:
steve
Message:

Adding old files as well as workinng files from anuga_1d

Location:
anuga_work/development/pipeflow
Files:
5 added
5 edited
2 copied

Legend:

Unmodified
Added
Removed
  • anuga_work/development/pipeflow/generic_domain.py

    r7788 r7823  
    1212
    1313
    14 class Generic_Domain:
     14class Generic_domain:
    1515
    1616    def __init__(self,
     
    587587            return self.vertices[elem_id,vertex]
    588588
    589     def get_area(self, elem_id):
     589    def get_area(self, elem_id=None):
    590590        """Return area of element id
    591591        """
    592592
    593         return self.areas[elem_id]
     593        if elem_id is None:
     594            return sum(self.areas)
     595        else:
     596            return self.areas[elem_id]
     597
    594598
    595599    def get_quantity(self, name, location='vertices', indices = None):
  • anuga_work/development/pipeflow/quantity.py

    r7784 r7823  
    1 """Class Quantity - Implements values at each 1d element
     1"""
     2Class Quantity - Implements values at each 1d element
    23
    34To create:
     
    1718import numpy
    1819
     20
     21
    1922class Quantity:
    2023
     
    2326        #Initialise Quantity using optional vertex values.
    2427       
    25         from generic_domain import Generic_Domain
     28        from generic_domain import Generic_domain
    2629
    2730        msg = 'First argument in Quantity.__init__ '
    2831        msg += 'must be of class Domain (or a subclass thereof)'
    29         assert isinstance(domain, Generic_Domain), msg
     32        assert isinstance(domain, Generic_domain), msg
    3033
    3134        if vertex_values is None:
     
    5255        self.centroid_values = numpy.zeros(N, numpy.float)
    5356        self.centroid_backup_values = numpy.zeros(N, numpy.float)
    54         #self.edge_values = zeros((N, 2), numpy.float)
     57        #self.edge_values = numpy.zeros((N, 2), numpy.float)
    5558        #edge values are values of the ends of each interval
    5659     
     
    5861        self.interpolate()
    5962
     63
    6064       
    6165        #Allocate space for boundary values
     
    8286        """
    8387        return self.centroid_values.shape[0]
    84    
     88
     89    def __neg__(self):
     90        """Negate all values in this quantity giving meaning to the
     91        expression -Q where Q is an instance of class Quantity
     92        """
     93
     94        Q = Quantity(self.domain)
     95        Q.set_values_from_numeric(-self.vertex_values)
     96        return Q
     97
     98
     99
     100    def __add__(self, other):
     101        """Add to self anything that could populate a quantity
     102
     103        E.g other can be a constant, an array, a function, another quantity
     104        (except for a filename or points, attributes (for now))
     105        - see set_values_from_numeric for details
     106        """
     107
     108        Q = Quantity(self.domain)
     109        Q.set_values_from_numeric(other)
     110
     111        result = Quantity(self.domain)
     112        result.set_values_from_numeric(self.vertex_values + Q.vertex_values)
     113        return result
     114
     115    def __radd__(self, other):
     116        """Handle cases like 7+Q, where Q is an instance of class Quantity
     117        """
     118
     119        return self + other
     120
     121    def __sub__(self, other):
     122        return self + -other            # Invoke self.__neg__()
     123
     124    def __mul__(self, other):
     125        """Multiply self with anything that could populate a quantity
     126
     127        E.g other can be a constant, an array, a function, another quantity
     128        (except for a filename or points, attributes (for now))
     129        - see set_values_from_numeric for details
     130        """
     131
     132        if isinstance(other, Quantity):
     133            Q = other
     134        else:
     135            Q = Quantity(self.domain)
     136            Q.set_values_from_numeric(other)
     137
     138        result = Quantity(self.domain)
     139
     140        # The product of vertex_values, edge_values and centroid_values
     141        # are calculated and assigned directly without using
     142        # set_values_from_numeric (which calls interpolate). Otherwise
     143        # centroid values wouldn't be products from q1 and q2
     144        result.vertex_values = self.vertex_values * Q.vertex_values
     145        result.centroid_values = self.centroid_values * Q.centroid_values
     146
     147        return result
     148
     149    def __rmul__(self, other):
     150        """Handle cases like 3*Q, where Q is an instance of class Quantity
     151        """
     152
     153        return self * other
     154
     155    def __div__(self, other):
     156        """Divide self with anything that could populate a quantity
     157
     158        E.g other can be a constant, an array, a function, another quantity
     159        (except for a filename or points, attributes (for now))
     160        - see set_values_from_numeric for details
     161
     162        Zero division is dealt with by adding an epsilon to the divisore
     163        FIXME (Ole): Replace this with native INF once we migrate to NumPy
     164        """
     165
     166        if isinstance(other, Quantity):
     167            Q = other
     168        else:
     169            Q = Quantity(self.domain)
     170            Q.set_values_from_numeric(other)
     171
     172        result = Quantity(self.domain)
     173
     174        # The quotient of vertex_values, edge_values and centroid_values
     175        # are calculated and assigned directly without using
     176        # set_values_from_numeric (which calls interpolate). Otherwise
     177        # centroid values wouldn't be quotient of q1 and q2
     178        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
     179        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
     180
     181        return result
     182
     183    def __rdiv__(self, other):
     184        """Handle cases like 3/Q, where Q is an instance of class Quantity
     185        """
     186
     187        return self / other
     188
     189    def __pow__(self, other):
     190        """Raise quantity to (numerical) power
     191
     192        As with __mul__ vertex values are processed entry by entry
     193        while centroid and edge values are re-interpolated.
     194
     195        Example using __pow__:
     196          Q = (Q1**2 + Q2**2)**0.5
     197        """
     198
     199        if isinstance(other, Quantity):
     200            Q = other
     201        else:
     202            Q = Quantity(self.domain)
     203            Q.set_values_from_numeric(other)
     204
     205        result = Quantity(self.domain)
     206
     207        # The power of vertex_values, edge_values and centroid_values
     208        # are calculated and assigned directly without using
     209        # set_values_from_numeric (which calls interpolate). Otherwise
     210        # centroid values wouldn't be correct
     211        result.vertex_values = self.vertex_values ** other
     212        result.centroid_values = self.centroid_values ** other
     213
     214        return result
     215
     216    def set_values_from_numeric(self, numeric):
     217       
     218        x = numpy.array([1.0,2.0])
     219        y = [1.0,2.0]
     220
     221        if type(numeric) == type(y):
     222            self.set_values_from_array(numeric)
     223        elif type(numeric) == type(x):
     224            self.set_values_from_array(numeric)
     225        elif callable(numeric):
     226            self.set_values_from_function(numeric)
     227        elif isinstance(numeric, Quantity):
     228            self.set_values_from_quantity(numeric)
     229        else:   # see if it's coercible to a float (float, int or long, etc)
     230            try:
     231                numeric = float(numeric)
     232            except ValueError:
     233                msg = ("Illegal type for variable 'numeric': %s" % type(numeric))
     234                raise Exception(msg)
     235            self.set_values_from_constant(numeric)
     236
     237    def set_values_from_constant(self,numeric):
     238
     239        self.vertex_values[:,:]   = numeric
     240        self.centroid_values[:,] = numeric
     241
     242
     243    def set_values_from_array(self,numeric):
     244
     245        self.vertex_values[:,:]   = numeric
     246        self.interpolate()
     247
     248
     249    def set_values_from_quantity(self,quantity):
     250
     251        self.vertex_values[:,:]   = quantity.vertex_values
     252        self.centroid_values[:,] = quantity.centroid_values
     253
     254    def set_values_from_function(self,function):
     255
     256        self.vertex_values[:,:]   = map(function, self.domain.vertices)
     257        self.centroid_values[:,] = map(function, self.domain.centroids)
     258
     259
    85260    def interpolate(self):
    86261        """
     
    96271            self.centroid_values[i] = (v0 + v1)/2.0
    97272
     273
     274
     275
    98276    def set_values(self, X, location='vertices'):
    99277        """Set values for quantity
    100278
    101         X: Compatible list, numpy array (see below), constant or function
     279        X: Compatible list, Numeric array (see below), constant or function
    102280        location: Where values are to be stored.
    103281                  Permissible options are: vertices, centroid
     
    105283
    106284        In case of location == 'centroid' the dimension values must
    107         be a list of a numpyal array of length N, N being the number
     285        be a list of a Numerical array of length N, N being the number
    108286        of elements in the mesh. Otherwise it must be of dimension Nx2
    109287
     
    119297        """
    120298
    121         if location not in ['vertices', 'centroids']:
    122             msg = 'Invalid location: %s, (possible choices vertices, centroids)' %location
    123             raise msg
     299        if location not in ['vertices', 'centroids', 'unique_vertices']:
     300            msg = 'Invalid location: %s, (possible choices vertices, centroids, unique_vertices)' %location
     301            raise Exception(msg)
    124302
    125303        if X is None:
    126304            msg = 'Given values are None'
    127             raise msg
     305            raise Exception(msg)
    128306
    129307        import types
     
    134312         
    135313        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
    136             if location == 'centroids':
    137                 self.centroid_values[:] = X
    138             else:
    139                 self.vertex_values[:] = X
     314           
     315            self.centroid_values[:,] = float(X)
     316            self.vertex_values[:,:] = float(X)
     317
     318        elif isinstance(X, Quantity):
     319            self.set_array_values(X.vertex_values, location)
    140320
    141321        else:
     
    143323            self.set_array_values(X, location)
    144324
    145         if location == 'vertices':
    146             #Intialise centroid and edge values
     325        if location == 'vertices' or location == 'unique_vertices':
     326            #Intialise centroid
    147327            self.interpolate()
     328
     329        if location == 'centroid':
     330            self.extrapolate_first_order()
    148331
    149332
     
    176359        """Set values for quantity
    177360
    178         values: numpy array
     361        values: Numeric array
    179362        location: Where values are to be stored.
    180                   Permissible options are: vertices, centroid, edges
     363                  Permissible options are: vertices, centroid, unique_vertices
    181364                  Default is "vertices"
    182365
    183366        In case of location == 'centroid' the dimension values must
    184         be a list of a numpyal array of length N, N being the number
    185         of elements in the mesh. Otherwise it must be of dimension Nx2
     367        be a list of a Numerical array of length N, N being the number
     368        of elements in the mesh. If location == 'unique_vertices' the
     369        dimension values must be a list or a Numeric array of length N+1.
     370        Otherwise it must be of dimension Nx2
    186371
    187372        The values will be stored in elements following their
     
    195380
    196381
    197         values = numpy.array(values, numpy.float)
     382        values = numpy.array(values).astype(numpy.float)
    198383
    199384        N = self.centroid_values.shape[0]
    200385
    201         msg = 'Number of values must match number of elements'
    202         assert values.shape[0] == N, msg
    203386
    204387        if location == 'centroids':
     388            msg = 'Number of values must match number of elements'
     389            assert values.shape[0] == N, msg
    205390            assert len(values.shape) == 1, 'Values array must be 1d'
    206             self.centroid_values = values
    207         #elif location == 'edges':
    208         #    assert len(values.shape) == 2, 'Values array must be 2d'
    209         #    msg = 'Array must be N x 2'
    210         #    self.edge_values = values
    211         else:
     391
     392            for i in range(N):
     393                self.centroid_values[i] = values[i]
     394
     395            self.vertex_values[:,0] = values
     396            self.vertex_values[:,1] = values
     397 
     398        if location == 'vertices':
     399            msg = 'Number of values must match number of elements'
     400            assert values.shape[0] == N, msg
    212401            assert len(values.shape) == 2, 'Values array must be 2d'
    213402            msg = 'Array must be N x 2'
    214403            assert values.shape[1] == 2, msg
    215404
    216             self.vertex_values = values
     405            self.vertex_values[:,:] = values[:,:]
     406
     407        if location == 'unique_vertices':
     408            msg = 'Number of values must match number of elements +1'
     409            assert values.shape[0] == N+1, msg
     410            assert len(values.shape) == 1, 'Values array must be 1d'
     411
     412            self.vertex_values[:,0] = values[0:-1]
     413            self.vertex_values[:,1] = values[1:N+1] 
     414
     415
     416           
    217417
    218418
     
    220420        """get values for quantity
    221421
    222         return X, Compatible list, numpy array (see below)
     422        return X, Compatible list, Numeric array (see below)
    223423        location: Where values are to be stored.
    224424                  Permissible options are: vertices, centroid
     
    226426
    227427        In case of location == 'centroids' the dimension values must
    228         be a list of a numpyal array of length N, N being the number
     428        be a list of a Numerical array of length N, N being the number
    229429        of elements. Otherwise it must be of dimension Nx3
    230430
     
    236436
    237437        """
    238         from numpy import take
     438
    239439
    240440        if location not in ['vertices', 'centroids', 'unique vertices']:
    241441            msg = 'Invalid location: %s' %location
    242             raise msg
    243 
    244         import types
     442            raise Exception(msg)
     443
     444        import types, numpy
    245445        assert type(indices) in [types.ListType, types.NoneType,
    246446                                 numpy.ArrayType],\
     
    262462                if cells is None:
    263463                    msg = 'Unique vertex not associated with cells'
    264                     raise msg
     464                    raise Exception(msg)
    265465
    266466                # Go through all cells, vertex pairs
     
    274474            if (indices ==  None):
    275475                indices = range(len(self))
    276             return numpy.take(self.vertex_values,indices)
     476            return take(self.vertex_values,indices)
    277477
    278478
     
    284484        """Return vertex values like an OBJ format
    285485
    286         The vertex values are returned as one sequence in the 1D numpy.float array A.
     486        The vertex values are returned as one sequence in the 1D float array A.
    287487        If requested the coordinates will be returned in 1D arrays X.
    288488
     
    308508
    309509        """
     510
     511
    310512
    311513
     
    399601        """
    400602
     603
     604
    401605        N = self.centroid_values.shape[0]
    402606
     
    407611        denominator = numpy.ones(N, numpy.float)-timestep*self.semi_implicit_update
    408612
    409         if numpy.sum(numpy.equal(denominator, 0.0)) > 0.0:
     613        if sum(numpy.equal(denominator, 0.0)) > 0.0:
    410614            msg = 'Zero division in semi implicit update. Call Stephen :-)'
    411             raise msg
     615            raise Exception(msg)
    412616        else:
    413617            #Update conserved_quantities from semi implicit updates
     
    420624        """
    421625
    422         #print 'compute_gradient'
     626
     627
    423628
    424629        N = self.centroid_values.shape[0]
     
    510715
    511716        #print 'compute_minmod_gradients'
     717        from numpy import sign
     718
    512719       
    513720        def xmin(a,b):
     
    516723        def xmic(t,a,b):
    517724            return xmin(t*xmin(a,b), 0.50*(a+b) )
     725
     726
    518727
    519728        N = self.centroid_values.shape[0]
     
    702911        """
    703912
     913
    704914        N = self.domain.number_of_elements
    705915        beta = self.domain.beta
     
    729939
    730940        #Deltas between vertex and centroid values
    731         dq = zeros(qv.shape, numpy.float)
     941        dq = numpy.zeros(qv.shape, numpy.float)
    732942        for i in range(2):
    733943            dq[:,i] = qv[:,i] - qc
     
    753963
    754964        from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
     965
    755966        limiter = self.domain.limiter
    756967        #print limiter
    757968       
    758969        #print 'limit_range'
    759         N  = self.domain.number_of_elements
     970        N = self.domain.number_of_elements
    760971        qc = self.centroid_values
    761972        qv = self.vertex_values
    762         C  = self.domain.centroids
    763         X  = self.domain.vertices
     973        C = self.domain.centroids
     974        X = self.domain.vertices
    764975        beta_p = numpy.zeros(N,numpy.float)
    765976        beta_m = numpy.zeros(N,numpy.float)
     
    777988                beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
    778989               
    779         dq = zeros(qv.shape, numpy.float)
     990        dq = numpy.zeros(qv.shape, numpy.float)
    780991        for i in range(2):
    781992            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
     
    8171028
    8181029        import sys
     1030
    8191031        from util import minmod, minmod_kurganov, maxmod, vanleer
    8201032
     
    8411053                # Check denominator not zero
    8421054                if (qc[k+1]-qc[k]) == 0.0:
    843                     beta_p[k] = numpy.float(sys.maxint)
    844                     beta_m[k] = numpy.float(sys.maxint)
     1055                    beta_p[k] = float(sys.maxint)
     1056                    beta_m[k] = float(sys.maxint)
    8451057                else:
    8461058                    #STEVE LIMIT
     
    8491061
    8501062        #Deltas between vertex and centroid values
    851         dq = zeros(qv.shape, numpy.float)
     1063        dq = numpy.zeros(qv.shape, numpy.float)
    8521064        for i in range(2):
    8531065            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
     
    8851097        #backup_centroid_values(self)
    8861098
    887         self.centroid_backup_values[:] = (self.centroid_values).astype('f')
     1099        self.centroid_backup_values[:,] = (self.centroid_values).astype('f')
    8881100
    8891101    def saxpy_centroid_values(self,a,b):
    8901102        # Call correct module function
    8911103        # (either from this module or C-extension)
    892         self.centroid_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
     1104        self.centroid_values[:,] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
    8931105       
    8941106
     
    9141126if __name__ == "__main__":
    9151127    #from domain import Domain
    916     from shallow_water_domain import Domain     
    917     from numpy import arange
     1128    from sww_domain import Domain
    9181129   
    9191130    points1 = [0.0, 1.0, 2.0, 3.0]
     
    9681179    raw_input('press return')
    9691180
    970     points2 = arange(10)
     1181    points2 = numpy.arange(10)
    9711182    D2 = Domain(points2)
    9721183
  • anuga_work/development/pipeflow/sww_comp_flux_ext.c

    r7783 r7823  
    1010
    1111
     12/* double max(double a, double b) { */
     13/*      double z; */
     14/*      z=(a>b)?a:b; */
     15/*      return z;} */
     16       
     17/* double min(double a, double b) { */
     18/*      double z; */
     19/*      z=(a<b)?a:b; */
     20/*      return z;} */
     21
     22
     23// Function to obtain speed from momentum and depth.
     24// This is used by flux functions
     25// Input parameters uh and h may be modified by this function.
     26double _compute_speed(double *uh,
     27                      double *h,
     28                      double epsilon,
     29                      double h0) {
     30 
     31  double u;
     32 
     33  if (*h < epsilon) {
     34    *h = 0.0;  //Could have been negative
     35     u = 0.0;
     36  } else {
     37    u = *uh/(*h + h0/ *h);   
     38  }
     39 
     40
     41  // Adjust momentum to be consistent with speed
     42  *uh = u * *h;
     43 
     44  return u;
     45}
     46
     47
     48
     49
    1250//Innermost flux function (using w=z+h)
    1351int _flux_function(double *q_left, double *q_right,
     52        double z_left, double z_right,
    1453                   double normals, double g, double epsilon, double h0,
    15                    double *edgeflux, double *max_speed) {
     54                double *edgeflux, double *max_speed) {
    1655               
    1756                int i;
    18                 double flux_left[2], flux_right[2];
    19                 double w_left, h_left, uh_left, z_left, u_left, soundspeed_left;
    20                 double w_right, h_right, uh_right, z_right, u_right, soundspeed_right;
    21                 double z, s_max, s_min, denom;
     57                double ql[2], qr[2], flux_left[2], flux_right[2];
     58                double z, w_left, h_left, uh_left, soundspeed_left, u_left;
     59                double w_right, h_right, uh_right, soundspeed_right, u_right;
     60                double s_max, s_min, denom;
    2261               
    2362                //printf("h0 = %f \n",h0);
    24                 w_left  = q_left[0];
    25                 uh_left = q_left[1]*normals;
    26                 h_left  = q_left[2];
    27                 z_left  = q_left[3];
    28                 u_left  = q_left[4]*normals;
    29        
    30                 w_right  = q_right[0];
    31                 uh_right = q_right[1]*normals;
    32                 h_right  = q_right[2];
    33                 z_right  = q_right[3];
    34                 u_right  = q_right[4]*normals;         
     63                ql[0] = q_left[0];
     64                ql[1] = q_left[1];
     65                ql[1] = ql[1]*normals;
     66       
     67                qr[0] = q_right[0];
     68                qr[1] = q_right[1];
     69                qr[1] = qr[1]*normals;
    3570         
    3671                z = (z_left+z_right)/2.0;                 
    3772                                   
     73                //w_left = ql[0];
     74                //h_left = w_left-z;
     75                //uh_left = ql[1];
     76               
     77
     78
     79                // Compute speeds in x-direction
     80                w_left = ql[0];         
     81                h_left = w_left-z;
     82                uh_left = ql[1];
     83
     84                u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
     85
     86                w_right = qr[0];
     87                h_right = w_right-z;
     88                uh_right = qr[1];
     89
     90                u_right = _compute_speed(&uh_right, &h_right, epsilon, h0);
     91 
    3892                soundspeed_left = sqrt(g*h_left);
    3993                soundspeed_right = sqrt(g*h_right);
     
    48102                // Flux formulas
    49103                flux_left[0] = u_left*h_left;
    50                 flux_left[1] = u_left*u_left*h_left + 0.5*g*h_left*h_left;
     104                flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
    51105
    52106                flux_right[0] = u_right*h_right;
    53                 flux_right[1] = u_right*u_right*h_right + 0.5*g*h_right*h_right;
     107                flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
    54108
    55109                // Flux computation
     
    60114                } else {
    61115                        edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
    62                         edgeflux[0] += s_max*s_min*(q_right[0]-q_left[0]);
     116                        edgeflux[0] += s_max*s_min*(qr[0]-ql[0]);
    63117                        edgeflux[0] /= denom;
    64118                        edgeflux[1] = s_max*flux_left[1] - s_min*flux_right[1];
    65                         edgeflux[1] += s_max*s_min*(q_right[1]-q_left[1]);
     119                        edgeflux[1] += s_max*s_min*(qr[1]-ql[1]);
    66120                        edgeflux[1] /= denom;
    67121                        edgeflux[1] *= normals;
    68122   
    69                         // Maximal wavespeed
    70                         *max_speed = max(fabs(s_max), fabs(s_min));
     123                // Maximal wavespeed
     124        *max_speed = max(fabs(s_max), fabs(s_min));
    71125                } 
    72                 return 0;               
    73 }
     126    return 0;           
     127        }
    74128               
    75129               
     
    90144                           double* xmom_edge_values,
    91145                           double* bed_edge_values,
    92                            double* height_edge_values,
    93                            double* velocity_edge_values,
    94146                           double* stage_boundary_values,
    95147                           double* xmom_boundary_values,
    96                            double* bed_boundary_values,
    97                            double* height_boundary_values,
    98                            double* velocity_boundary_values,
    99148                           double* stage_explicit_update,
    100149                           double* xmom_explicit_update,
     
    102151                           double* max_speed_array) {
    103152               
    104                 double flux[2], ql[5], qr[5], edgeflux[2];
    105                 double max_speed, normal;
     153                double flux[2], ql[2], qr[2], edgeflux[2];
     154                double zl, zr, max_speed, normal;
    106155                int k, i, ki, n, m, nm=0;
    107156               
     
    116165                                ql[0] = stage_edge_values[ki];
    117166                                ql[1] = xmom_edge_values[ki];
    118                                 ql[2] = bed_edge_values[ki];
    119                                 ql[3] = height_edge_values[ki];
    120                                 ql[4] = velocity_edge_values[ki];
     167                                zl = bed_edge_values[ki];
    121168                               
    122169                                n = neighbours[ki];
     
    125172                                        qr[0] = stage_boundary_values[m];
    126173                                        qr[1] = xmom_boundary_values[m];
    127                                         qr[2] = bed_edge_values[m];
    128                                         qr[3] = height_edge_values[m];
    129                                         qr[4] = velocity_edge_values[m];
     174                                        zr = zl;
    130175                                } else {
    131176                                        m = neighbour_vertices[ki];
     
    133178                                        qr[0] = stage_edge_values[nm];
    134179                                        qr[1] = xmom_edge_values[nm];
    135                                         qr[2] = bed_edge_values[nm];
    136                                         qr[3] = height_edge_values[nm];
    137                                         qr[4] = velocity_edge_values[nm];
     180                                        zr = bed_edge_values[nm];                               
    138181                                }
    139182                               
    140183                                normal = normals[ki];
    141                                 _flux_function(ql, qr, normal, g, epsilon, h0, edgeflux, &max_speed);
     184                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed);
    142185                                flux[0] -= edgeflux[0];
    143186                                flux[1] -= edgeflux[1];
     
    152195                                    }
    153196                                }
    154                         } // End edge i (and neighbour n)
     197            } // End edge i (and neighbour n)
    155198                        flux[0] /= areas[k];
    156199                        stage_explicit_update[k] = flux[0];
     
    163206                return timestep;       
    164207        }
     208       
     209       
     210       
     211       
     212       
     213
     214
    165215
    166216//=========================================================================
     
    168218//=========================================================================
    169219PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
    170 
    171     PyObject
     220 
     221    PyArrayObject
     222        *neighbours,
     223        *neighbour_vertices,
     224        *normals,
     225        *areas,
     226        *stage_edge_values,
     227        *xmom_edge_values,
     228        *bed_edge_values,
     229        *stage_boundary_values,
     230        *xmom_boundary_values,
     231        *stage_explicit_update,
     232        *xmom_explicit_update,
     233        *max_speed_array;
     234   
     235  double timestep, epsilon, g, h0, cfl;
     236  int number_of_elements;
     237   
     238  // Convert Python arguments to C
     239  if (!PyArg_ParseTuple(args, "dddddOOOOOOOOOOOiO",
     240                        &cfl,
     241                        &timestep,
     242                        &epsilon,
     243                        &g,
     244                        &h0,
     245                        &neighbours,
     246                        &neighbour_vertices,
     247                        &normals,
     248                        &areas,
     249                        &stage_edge_values,
     250                        &xmom_edge_values,
     251                        &bed_edge_values,
     252                        &stage_boundary_values,
     253                        &xmom_boundary_values,
     254                        &stage_explicit_update,
     255                        &xmom_explicit_update,
     256                        &number_of_elements,
     257                        &max_speed_array)) {
     258    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
     259    return NULL;
     260  }
     261
     262 
     263  // Call underlying flux computation routine and update
     264  // the explicit update arrays
     265  timestep = _compute_fluxes_ext(cfl,
     266                                 timestep,
     267                                 epsilon,
     268                                 g,
     269                                 h0,
     270                                 (long*) neighbours -> data,
     271                                 (long*) neighbour_vertices -> data,
     272                                 (double*) normals -> data,
     273                                 (double*) areas -> data,
     274                                 (double*) stage_edge_values -> data,
     275                                 (double*) xmom_edge_values -> data,
     276                                 (double*) bed_edge_values -> data,
     277                                 (double*) stage_boundary_values -> data,
     278                                 (double*) xmom_boundary_values -> data,
     279                                 (double*) stage_explicit_update -> data,
     280                                 (double*) xmom_explicit_update -> data,
     281                                 number_of_elements,
     282                                 (double*) max_speed_array -> data);
     283
     284  // Return updated flux timestep
     285  return Py_BuildValue("d", timestep);
     286}
     287
     288
     289//-----------------------
     290PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) {
     291 
     292   PyObject
    172293        *domain,
    173         *stage,
    174         *xmom,
    175         *bed,
    176         *height,
    177         *velocity;
    178 
    179     PyArrayObject
    180         *neighbours,
     294        *stage,
     295        *xmom,
     296        *bed;
     297
     298    PyArrayObject
     299        *neighbours,
    181300        *neighbour_vertices,
    182         *normals,
     301        *normals, 
    183302        *areas,
    184303        *stage_vertex_values,
    185304        *xmom_vertex_values,
    186305        *bed_vertex_values,
    187         *height_vertex_values,
    188         *velocity_vertex_values,
    189306        *stage_boundary_values,
    190307        *xmom_boundary_values,
    191         *bed_boundary_values,
    192         *height_boundary_values,
    193         *velocity_boundary_values,
    194308        *stage_explicit_update,
    195309        *xmom_explicit_update,
    196310        *max_speed_array;
    197 
     311   
    198312  double timestep, epsilon, g, h0, cfl;
    199313  int number_of_elements;
    200314
    201 
     315   
    202316  // Convert Python arguments to C
    203   if (!PyArg_ParseTuple(args, "dOOOOOO",
     317  if (!PyArg_ParseTuple(args, "dOOOO",
    204318                        &timestep,
    205319                        &domain,
    206320                        &stage,
    207321                        &xmom,
    208                         &bed,
    209                         &height,
    210                         &velocity)) {
     322                        &bed)) {
    211323      PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
    212324      return NULL;
     
    217329    g                 = get_python_double(domain,"g");
    218330    h0                = get_python_double(domain,"h0");
    219     cfl               = get_python_double(domain,"cfl");
    220 
    221 
     331    cfl               = get_python_double(domain,"CFL");
     332 
     333   
    222334    neighbours        = get_consecutive_array(domain, "neighbours");
    223     neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
     335    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices"); 
    224336    normals           = get_consecutive_array(domain, "normals");
    225     areas             = get_consecutive_array(domain, "areas");
     337    areas             = get_consecutive_array(domain, "areas");   
    226338    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
    227 
    228     stage_vertex_values      = get_consecutive_array(stage,    "vertex_values");
    229     xmom_vertex_values       = get_consecutive_array(xmom,     "vertex_values");
    230     bed_vertex_values        = get_consecutive_array(bed,      "vertex_values");
    231     height_vertex_values     = get_consecutive_array(height,   "vertex_values");
    232     velocity_vertex_values   = get_consecutive_array(velocity, "vertex_values");
    233 
    234     stage_boundary_values     = get_consecutive_array(stage,     "boundary_values");
    235     xmom_boundary_values      = get_consecutive_array(xmom,      "boundary_values");
    236     bed_boundary_values       = get_consecutive_array(bed,       "boundary_values");
    237     height_boundary_values    = get_consecutive_array(height,    "boundary_values");
    238     velocity_boundary_values  = get_consecutive_array(velocity,  "boundary_values");
    239 
    240 
    241     stage_explicit_update = get_consecutive_array(stage, "explicit_update");
    242     xmom_explicit_update  = get_consecutive_array(xmom,  "explicit_update");
     339   
     340    stage_vertex_values = get_consecutive_array(stage, "vertex_values");   
     341    xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");   
     342    bed_vertex_values   = get_consecutive_array(bed, "vertex_values");   
     343
     344    stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
     345    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
     346
     347
     348    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
     349    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
     350
     351
    243352
    244353    number_of_elements = stage_vertex_values -> dimensions[0];
    245354
    246     // Call underlying flux computation routine and update
    247     // the explicit update arrays
     355
     356 
     357    // Call underlying flux computation routine and update
     358    // the explicit update arrays
    248359    timestep = _compute_fluxes_ext(
    249360        cfl,
     
    259370        (double*) xmom_vertex_values -> data,
    260371        (double*) bed_vertex_values -> data,
    261         (double*) height_vertex_values -> data,
    262         (double*) velocity_vertex_values -> data,
    263372        (double*) stage_boundary_values -> data,
    264373        (double*) xmom_boundary_values -> data,
    265         (double*) bed_boundary_values -> data,
    266         (double*) height_boundary_values -> data,
    267         (double*) velocity_boundary_values -> data,
    268374        (double*) stage_explicit_update -> data,
    269375        (double*) xmom_explicit_update -> data,
     
    279385  Py_DECREF(xmom_vertex_values);
    280386  Py_DECREF(bed_vertex_values);
    281   Py_DECREF(height_vertex_values);
    282   Py_DECREF(velocity_vertex_values);
    283387  Py_DECREF(stage_boundary_values);
    284388  Py_DECREF(xmom_boundary_values);
    285   Py_DECREF(bed_boundary_values);
    286   Py_DECREF(height_boundary_values);
    287   Py_DECREF(velocity_boundary_values);
    288389  Py_DECREF(stage_explicit_update);
    289390  Py_DECREF(xmom_explicit_update);
     
    291392
    292393
     394
     395
    293396  // Return updated flux timestep
    294397  return Py_BuildValue("d", timestep);
    295398}
    296399
     400 
    297401
    298402
     
    303407static struct PyMethodDef MethodTable[] = {
    304408  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
     409  {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
    305410  {NULL, NULL}
    306411};
  • anuga_work/development/pipeflow/sww_domain.py

    r7781 r7823  
    5050
    5151#Shallow water domain
    52 class Domain(Generic_Domain):
     52class Domain(Generic_domain):
    5353
    5454    def __init__(self, coordinates, boundary = None, tagged_elements = None):
     
    5757        evolved_quantities   = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
    5858        other_quantities     = ['friction']
    59         Generic_Domain.__init__(self,
     59        Generic_domain.__init__(self,
    6060                                coordinates          = coordinates,
    6161                                boundary             = boundary,
     
    143143
    144144
    145         timestep = numpy.float(sys.maxint)
     145        timestep = float(sys.maxint)
    146146
    147147        stage    = self.quantities['stage']
     
    152152
    153153
    154         from sww_comp_flux_ext import compute_fluxes_ext
    155 
    156         self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity)
     154        from sww_comp_flux_ext import compute_fluxes_ext_short
     155
     156        #self.flux_timestep = compute_fluxes_ext(timestep,self,stage,xmom,bed,height,velocity)
     157
     158   
     159        self.flux_timestep = compute_fluxes_ext_short(timestep,self,stage,xmom,bed)
     160
    157161
    158162
     
    168172
    169173        #Call basic machinery from parent class
    170         for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
     174        for t in Generic_domain.evolve(self, yieldstep, finaltime,duration,
    171175                                       skip_initial_step):
    172176
     
    266270##      else:
    267271
    268     u_C[:]  = uh_C/(h_C + h0/h_C)
     272    u_C[:,]  = uh_C/(h_C + h0/h_C)
    269273       
    270274    for name in [ 'velocity', 'stage' ]:
     
    286290    xmom_V  = domain.quantities['xmomentum'].vertex_values
    287291
    288     h_V[:]    = stage_V - bed_V
     292    h_V[:,:]    = stage_V - bed_V
    289293    for i in range(len(h_C)):
    290294        for j in range(2):
     
    297301                stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
    298302               
    299     xmom_V[:] = u_V * h_V
     303    xmom_V[:,:] = u_V * h_V
    300304   
    301305    return
     
    314318
    315319    #Shortcuts
    316     wc = domain.quantities['stage'].centroid_values
    317     zc = domain.quantities['elevation'].centroid_values
     320    wc    = domain.quantities['stage'].centroid_values
     321    zc    = domain.quantities['elevation'].centroid_values
    318322    xmomc = domain.quantities['xmomentum'].centroid_values
    319323    hc = wc - zc  #Water depths at centroids
  • anuga_work/development/pipeflow/sww_python.py

    r7789 r7823  
    5858    for k in range(N):
    5959
    60         flux[:] = 0.  #Reset work numpy.array
     60        flux[:,] = 0.  #Reset work numpy.array
    6161        #for i in range(3):
    6262        for i in range(2):
Note: See TracChangeset for help on using the changeset viewer.