Changeset 242


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

Subclassed class Quantity into Conserved_quantity

Location:
inundation/ga/storm_surge/pyvolution
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/config.py

    r240 r242  
    6262
    6363use_psyco = True  #Use psyco optimisations
    64 #use_psyco = False  #Do not use psyco optimisations
     64use_psyco = False  #Do not use psyco optimisations
    6565
    6666
  • inundation/ga/storm_surge/pyvolution/domain.py

    r240 r242  
    1919
    2020        from Numeric import zeros, Float
    21         from quantity import Quantity
     21        from quantity import Quantity, Conserved_quantity
    2222
    2323        #List of quantity names entering
     
    3434        #Build dictionary of Quantity instances keyed by quantity names
    3535        self.quantities = {}
    36         for name in self.conserved_quantities + other_quantities:
     36        for name in self.conserved_quantities:
     37            self.quantities[name] = Conserved_quantity(self)
     38        for name in other_quantities:
    3739            self.quantities[name] = Quantity(self)
    3840
  • inundation/ga/storm_surge/pyvolution/quantity.py

    r234 r242  
    1515"""
    1616
    17 #FIXME: Make Conserved_quantity a subclass of Quantity
    18 
    19 from mesh import Mesh
    2017
    2118class Quantity:
     
    2320    def __init__(self, domain, vertex_values=None):
    2421
     22        from mesh import Mesh
    2523        from Numeric import array, zeros, Float
    2624
     
    5351        self.edge_values = zeros((N, 3), Float)       
    5452
    55         #Allocate space for boundary values
    56         L = len(domain.boundary)
    57         self.boundary_values = zeros(L, Float)
    58 
    59         #Allocate space for updates of conserved quantities by
    60         #flux calculations and forcing functions
    61         self.explicit_update = zeros(N, Float )
    62         self.semi_implicit_update = zeros(N, Float )
    63                
    6453        #Intialise centroid and edge_values
    6554        self.interpolate()
     
    7160        """
    7261
    73         #FIXME: Write using vector operations (or in C)
    7462        N = self.vertex_values.shape[0]       
    7563        for i in range(N):
     
    7967           
    8068            self.centroid_values[i] = (v0 + v1 + v2)/3
    81            
    82             #FIXME: This is a duplicate from 'interpolate_from_vertices_to_edges'
    83             self.edge_values[i, 0] = 0.5*(v1 + v2)
    84             self.edge_values[i, 1] = 0.5*(v0 + v2)           
    85             self.edge_values[i, 2] = 0.5*(v0 + v1)
    86 
     69
     70        self.interpolate_from_vertices_to_edges()
     71
     72
     73    def interpolate_from_vertices_to_edges(self):
     74        for k in range(self.vertex_values.shape[0]):
     75            q0 = self.vertex_values[k, 0]
     76            q1 = self.vertex_values[k, 1]
     77            q2 = self.vertex_values[k, 2]
     78           
     79            self.edge_values[k, 0] = 0.5*(q1+q2)
     80            self.edge_values[k, 1] = 0.5*(q0+q2) 
     81            self.edge_values[k, 2] = 0.5*(q0+q1)
     82
     83           
     84           
     85    def set_values(self, X, location='vertices'):
     86        """Set values for quantity
     87
     88        X: Compatible list, Numeric array (see below), constant or function
     89        location: Where values are to be stored.
     90                  Permissible options are: vertices, edges, centroid
     91                  Default is "vertices"
     92
     93        In case of location == 'centroid' the dimension values must
     94        be a list of a Numerical array of length N, N being the number
     95        of elements in the mesh. Otherwise it must be of dimension Nx3
     96
     97        The values will be stored in elements following their
     98        internal ordering.
     99
     100        If values are described a function, it will be evaluated at specified points
     101
     102        If selected location is vertices, values for centroid and edges
     103        will be assigned interpolated values.
     104        In any other case, only values for the specified locations
     105        will be assigned and the others will be left undefined.
     106        """
     107
     108        if location not in ['vertices', 'centroids', 'edges']:
     109            msg = 'Invalid location: %s' %location
     110            raise msg
     111
     112        if X is None:
     113            msg = 'Given values are None'
     114            raise msg           
     115       
     116        import types
     117       
     118        if callable(X):
     119            #Use function specific method
     120            self.set_function_values(X, location)           
     121        elif type(X) in [types.FloatType, types.IntType, types.LongType]:
     122            if location == 'centroids':
     123                self.centroid_values[:] = X
     124            elif location == 'edges':
     125                self.edge_values[:] = X
     126            else:
     127                self.vertex_values[:] = X               
     128
     129        else:
     130            #Use array specific method
     131            self.set_array_values(X, location)
     132
     133        if location == 'vertices':
     134            #Intialise centroid and edge_values
     135            self.interpolate()
     136           
     137
     138
     139    def set_function_values(self, f, location='vertices'):
     140        """Set values for quantity using specified function
     141
     142        f: x, y -> z Function where x, y and z are arrays
     143        location: Where values are to be stored.
     144                  Permissible options are: vertices, edges, centroid
     145                  Default is "vertices"
     146        """
     147
     148        if location == 'centroids':
     149            P = self.domain.centroids
     150            self.set_values(f(P[:,0], P[:,1]), location)
     151        elif location == 'edges':
     152            raise 'Not implemented: %s' %location
     153        else:
     154            #Vertices
     155            P = self.domain.get_vertex_coordinates()
     156            for i in range(3):
     157                self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1])
     158
     159           
     160    def set_array_values(self, values, location='vertices'):
     161        """Set values for quantity
     162
     163        values: Numeric array
     164        location: Where values are to be stored.
     165                  Permissible options are: vertices, edges, centroid
     166                  Default is "vertices"
     167
     168        In case of location == 'centroid' the dimension values must
     169        be a list of a Numerical array of length N, N being the number
     170        of elements in the mesh. Otherwise it must be of dimension Nx3
     171
     172        The values will be stored in elements following their
     173        internal ordering.
     174
     175        If selected location is vertices, values for centroid and edges
     176        will be assigned interpolated values.
     177        In any other case, only values for the specified locations
     178        will be assigned and the others will be left undefined.
     179        """
     180
     181        from Numeric import array, Float
     182
     183        values = array(values).astype(Float)
     184
     185        N = self.centroid_values.shape[0]
     186       
     187        msg = 'Number of values must match number of elements'
     188        assert values.shape[0] == N, msg
     189       
     190        if location == 'centroids':
     191            assert len(values.shape) == 1, 'Values array must be 1d'
     192            self.centroid_values = values                       
     193        elif location == 'edges':
     194            assert len(values.shape) == 2, 'Values array must be 2d'
     195            msg = 'Array must be N x 3'           
     196            assert values.shape[1] == 3, msg
     197           
     198            self.edge_values = values
     199        else:
     200            assert len(values.shape) == 2, 'Values array must be 2d'
     201            msg = 'Array must be N x 3'           
     202            assert values.shape[1] == 3, msg
     203           
     204            self.vertex_values = values
     205
     206
     207
     208class Conserved_quantity(Quantity):
     209    """Class conserved quantity adds to Quantity:
     210
     211    boundary values, storage and method for updating, and
     212    methods for extrapolation from centropid to vertices inluding
     213    gradients and limiters
     214    """
     215
     216    def __init__(self, domain, vertex_values=None):
     217        Quantity.__init__(self, domain, vertex_values)
     218       
     219        from Numeric import zeros, Float
     220
     221        #Allocate space for boundary values
     222        L = len(domain.boundary)
     223        self.boundary_values = zeros(L, Float)
     224
     225        #Allocate space for updates of conserved quantities by
     226        #flux calculations and forcing functions
     227
     228        N = domain.number_of_elements
     229        self.explicit_update = zeros(N, Float )
     230        self.semi_implicit_update = zeros(N, Float )
     231               
    87232
    88233    def update(self, timestep):
     
    134279                #We have zero neighbouring volumes -
    135280                #Fall back to first order scheme
    136 
    137281                pass
    138282           
     
    162306                    a[k] = (y1*q0 - y0*q1)/det
    163307                    b[k] = (x0*q1 - x1*q0)/det
    164 
    165308
    166309            else:
     
    288431
    289432           
    290 
    291     def interpolate_from_vertices_to_edges(self):
    292         #FIXME: Write using vector operations (or in C)
    293         for k in range(self.vertex_values.shape[0]):
    294             q0 = self.vertex_values[k, 0]
    295             q1 = self.vertex_values[k, 1]
    296             q2 = self.vertex_values[k, 2]
    297            
    298             self.edge_values[k, 0] = 0.5*(q1+q2)
    299             self.edge_values[k, 1] = 0.5*(q0+q2) 
    300             self.edge_values[k, 2] = 0.5*(q0+q1)
    301 
    302            
    303            
    304     def set_values(self, X, location='vertices'):
    305         """Set values for quantity
    306 
    307         X: Compatible list, Numeric array (see below), constant or function
    308         location: Where values are to be stored.
    309                   Permissible options are: vertices, edges, centroid
    310                   Default is "vertices"
    311 
    312         In case of location == 'centroid' the dimension values must
    313         be a list of a Numerical array of length N, N being the number
    314         of elements in the mesh. Otherwise it must be of dimension Nx3
    315 
    316         The values will be stored in elements following their
    317         internal ordering.
    318 
    319         If values are described a function, it will be evaluated at specified points
    320 
    321         If selected location is vertices, values for centroid and edges
    322         will be assigned interpolated values.
    323         In any other case, only values for the specified locations
    324         will be assigned and the others will be left undefined.
    325         """
    326 
    327 
    328         #FIXME: Should do location error check here only
    329        
    330         if location not in ['vertices', 'centroids', 'edges']:
    331             msg = 'Invalid location: %s' %location
    332             raise msg
    333 
    334         if X is None:
    335             msg = 'Given values are None'
    336             raise msg           
    337        
    338         import types
    339        
    340         if callable(X):
    341             #Use function specific method
    342             self.set_function_values(X, location)           
    343         elif type(X) in [types.FloatType, types.IntType, types.LongType]:
    344             if location == 'centroids':
    345                 self.centroid_values[:] = X
    346             elif location == 'edges':
    347                 self.edge_values[:] = X
    348             else:
    349                 self.vertex_values[:] = X               
    350 
    351         else:
    352             #Use array specific method
    353             self.set_array_values(X, location)
    354 
    355         if location == 'vertices':
    356             #Intialise centroid and edge_values
    357             self.interpolate()
    358            
    359 
    360 
    361     def set_function_values(self, f, location='vertices'):
    362         """Set values for quantity using specified function
    363 
    364         f: x, y -> z Function where x, y and z are arrays
    365         location: Where values are to be stored.
    366                   Permissible options are: vertices, edges, centroid
    367                   Default is "vertices"
    368         """
    369 
    370         if location == 'centroids':
    371             P = self.domain.centroids
    372             self.set_values(f(P[:,0], P[:,1]), location)
    373         elif location == 'edges':
    374             raise 'Not yet implemented: %s' %location
    375         else:
    376             #Vertices
    377             P = self.domain.get_vertex_coordinates()
    378             for i in range(3):
    379                 self.vertex_values[:,i] = f(P[:,2*i], P[:,2*i+1])
    380 
    381            
    382     def set_array_values(self, values, location='vertices'):
    383         """Set values for quantity
    384 
    385         values: Numeric array
    386         location: Where values are to be stored.
    387                   Permissible options are: vertices, edges, centroid
    388                   Default is "vertices"
    389 
    390         In case of location == 'centroid' the dimension values must
    391         be a list of a Numerical array of length N, N being the number
    392         of elements in the mesh. Otherwise it must be of dimension Nx3
    393 
    394         The values will be stored in elements following their
    395         internal ordering.
    396 
    397         If selected location is vertices, values for centroid and edges
    398         will be assigned interpolated values.
    399         In any other case, only values for the specified locations
    400         will be assigned and the others will be left undefined.
    401         """
    402 
    403         from Numeric import array, Float
    404 
    405         values = array(values).astype(Float)
    406 
    407         N = self.centroid_values.shape[0]
    408        
    409         msg = 'Number of values must match number of elements'
    410         assert values.shape[0] == N, msg
    411        
    412         if location == 'centroids':
    413             assert len(values.shape) == 1, 'Values array must be 1d'
    414             self.centroid_values = values                       
    415         elif location == 'edges':
    416             assert len(values.shape) == 2, 'Values array must be 2d'
    417             msg = 'Array must be N x 3'           
    418             assert values.shape[1] == 3, msg
    419            
    420             self.edge_values = values
    421         else:
    422             assert len(values.shape) == 2, 'Values array must be 2d'
    423             msg = 'Array must be N x 3'           
    424             assert values.shape[1] == 3, msg
    425            
    426             self.vertex_values = values
    427 
  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r240 r242  
    352352        Q.interpolate_from_vertices_to_edges()           
    353353
     354
    354355def protect_against_infinitesimal_heights(domain):
    355356    """Protect against infinitesimal heights and high speeds   
     
    397398    #Update
    398399    for k in range(domain.number_of_elements):
    399 
    400 
    401400        if domain.newstyle:
    402401            if hc[k] < domain.minimum_allowed_height:       
  • inundation/ga/storm_surge/pyvolution/test_quantity.py

    r212 r242  
    8787
    8888    def test_boundary_allocation(self):
    89         quantity = Quantity(self.mesh4,
     89        quantity = Conserved_quantity(self.mesh4,
    9090                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    9191
     
    189189
    190190    def test_gradient(self):
    191         quantity = Quantity(self.mesh4)
     191        quantity = Conserved_quantity(self.mesh4)
    192192
    193193        #Set up for a gradient of (3,0) at mid triangle
     
    218218
    219219    def test_second_order_extrapolation2(self):
    220         quantity = Quantity(self.mesh4)
     220        quantity = Conserved_quantity(self.mesh4)       
    221221
    222222        #Set up for a gradient of (3,1), f(x) = 3x+y
     
    327327
    328328    def test_first_order_extrapolator(self):
    329         quantity = Quantity(self.mesh4)
     329        quantity = Conserved_quantity(self.mesh4)               
    330330
    331331        #Test centroids
     
    342342
    343343    def test_second_order_extrapolator(self):
    344         quantity = Quantity(self.mesh4)
     344        quantity = Conserved_quantity(self.mesh4)                       
    345345
    346346        #Set up for a gradient of (3,0) at mid triangle
     
    375375       
    376376
    377     def test_limiter(self):       
    378         quantity = Quantity(self.mesh4)
     377    def test_limiter(self):
     378        quantity = Conserved_quantity(self.mesh4)
    379379
    380380        #Create a deliberate overshoot (e.g. from gradient computation)
     
    406406
    407407    def test_distribute_first_order(self):
    408         quantity = Quantity(self.mesh4)
     408        quantity = Conserved_quantity(self.mesh4)       
    409409
    410410        #Test centroids
     
    427427
    428428    def test_update_explicit(self):
    429         quantity = Quantity(self.mesh4)
     429        quantity = Conserved_quantity(self.mesh4)
    430430
    431431        #Test centroids
     
    443443
    444444    def test_update_semi_implicit(self):
    445         quantity = Quantity(self.mesh4)
     445        quantity = Conserved_quantity(self.mesh4)
    446446
    447447        #Test centroids
     
    459459
    460460    def test_both_updates(self):
    461         quantity = Quantity(self.mesh4)
     461        quantity = Conserved_quantity(self.mesh4)
    462462
    463463        #Test centroids
Note: See TracChangeset for help on using the changeset viewer.