Changeset 7562 for anuga_core/source


Ignore:
Timestamp:
Nov 19, 2009, 5:23:52 PM (14 years ago)
Author:
steve
Message:

Updating the balanced and parallel code

Location:
anuga_core/source
Files:
1 added
14 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r7519 r7562  
    11381138        # Input checks
    11391139        if quantities is None:
    1140             quantities = self.conserved_quantities
     1140            quantities = self.evolved_quantities
    11411141        elif type(quantities) == types.StringType:
    11421142            quantities = [quantities] #Turn it into a list
     
    17521752    ##
    17531753    # @brief Mapping between conserved quantites and evolved quantities
    1754     # @param q_cons array of conserved quantity values
    1755     # @param q_evolved array of current evolved quantity values
    1756     def  conserved_to_evolved(self, q_cons, q_evolved):
     1754    # @param Input: q_cons array of conserved quantity values
     1755    # @param Input: q_evol array of current evolved quantity values
     1756    # @note  Output: Updated q_evol array
     1757    def  conserved_values_to_evolved_values(self, q_cons, q_evol):
    17571758        """Needs to be overridden by Domain subclass
    17581759        """
    17591760
    1760         if len(q_cons) == len(q_evolved):
    1761             q_evolved[:] = q_cons
     1761        if len(q_cons) == len(q_evol):
     1762            q_evol[:] = q_cons
    17621763        else:
    1763             msg = 'Method conserved_to_evolved must be overridden by Domain subclass'
     1764            msg = 'Method conserved_values_to_evolved_values must be overridden by Domain subclass'
    17641765            raise Exception, msg
     1766
     1767        return q_evol
    17651768   
    17661769    ##
     
    17811784                log.critical('WARNING: Ignored boundary segment (None)')
    17821785            else:
    1783                 q_cons = B.evaluate(vol_id, edge_id)
    1784 
    1785                 if len(q_cons) == len(self.evolved_quantities):
     1786                q_bdry = B.evaluate(vol_id, edge_id)
     1787
     1788                if len(q_bdry) == len(self.evolved_quantities):
    17861789                    # conserved and evolved quantities are the same
    1787                     q_evol = q_cons
    1788                 elif len(q_cons) == len(self.conserved_quantities):
     1790                    q_evol = q_bdry
     1791                elif len(q_bdry) == len(self.conserved_quantities):
    17891792                    # boundary just returns conserved quantities
    17901793                    # Need to calculate all the evolved quantities
     
    17931796                    q_evol = self.get_evolved_quantities(vol_id, edge = edge_id)
    17941797
    1795                     self.conserved_to_evolved(q_cons, q_evol)
     1798                    q_evol = self.conserved_values_to_evolved_values(q_bdry, q_evol)
    17961799                else:
    17971800                    msg = 'Boundary must return array of either conserved or evolved quantities'
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

    r7519 r7562  
    606606            pass
    607607        else:
    608             msg = 'Should have caught the lack of conserved_to_evolved member function'
     608            msg = 'Should have caught the lack of conserved_values_to_evolved_values member function'
    609609            raise Exception, msg
    610610
    611611
    612         def  conserved_to_evolved(q_cons, q_evol):
     612        def  conserved_values_to_evolved_values(q_cons, q_evol):
    613613
    614614            q_evol[0:3] = q_cons
     
    616616            q_evol[4] = q_cons[2]/q_cons[0]
    617617
    618         domain.conserved_to_evolved = conserved_to_evolved
     618            return q_evol
     619
     620        domain.conserved_values_to_evolved_values = conserved_values_to_evolved_values
    619621
    620622        domain.update_boundary()
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_util.py

    r7342 r7562  
    140140        domain1.store = True
    141141        domain1.set_datadir('.')
    142         domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
     142        sww_file = 'spatio_temporal_boundary_source_%d' %(id(self))
     143        domain1.set_name(sww_file)
    143144
    144145        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
  • anuga_core/source/anuga/culvert_flows/test_culvert_class.py

    r7497 r7562  
    124124       
    125125
     126    os.remove('Test_culvert_shallow.sww')
    126127
    127128class Test_Culvert(unittest.TestCase):
     
    243244        assert min_delta_w < 0
    244245        assert max_delta_w > 10       
    245        
     246
     247
     248        os.remove('Test_culvert.sww')
    246249
    247250    def test_that_culvert_dry_bed_rating_does_not_produce_flow(self):
     
    346349
    347350   
    348 
     351        os.remove('Test_culvert_dry.sww')
    349352       
    350353    def test_that_culvert_flows_conserves_volume(self):
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r7520 r7562  
    644644            fid.close()
    645645
     646
     647##
     648# @brief Class to open an sww file so that domain can be populated with quantity values
     649class Read_sww:
     650
     651    def __init__(self, source):
     652        """The source parameter is assumed to be a NetCDF sww file.
     653        """
     654
     655        self.source = source
     656
     657        self.frame_number = 0
     658
     659        fin = NetCDFFile(self.source, 'r')
     660
     661        self.time = num.array(fin.variables['time'], num.float)
     662        self.last_frame_number = self.time.shape[0] - 1
     663
     664        self.frames = num.arange(self.last_frame_number+1)
     665
     666        fin.close()
     667       
     668        self.read_mesh()
     669
     670        self.quantities = {}
     671
     672        self.read_quantities()
     673
     674
     675    def read_mesh(self):
     676        fin = NetCDFFile(self.source, 'r')
     677
     678        self.vertices = num.array(fin.variables['volumes'], num.int)
     679       
     680        self.x = x = num.array(fin.variables['x'], num.float)
     681        self.y = y = num.array(fin.variables['y'], num.float)
     682
     683        assert len(self.x) == len(self.y)
     684       
     685        self.xmin = num.min(x)
     686        self.xmax = num.max(x)
     687        self.ymin = num.min(y)
     688        self.ymax = num.max(y)
     689
     690
     691
     692        fin.close()
     693       
     694    def read_quantities(self, frame_number=0):
     695
     696        assert frame_number >= 0 and frame_number <= self.last_frame_number
     697
     698        self.frame_number = frame_number
     699
     700        M = len(self.x)/3
     701       
     702        fin = NetCDFFile(self.source, 'r')
     703       
     704        for q in filter(lambda n:n != 'x' and n != 'y' and n != 'z' and n != 'time' and n != 'volumes' and \
     705                        '_range' not in n, \
     706                        fin.variables.keys()):
     707            if len(fin.variables[q].shape) == 1: # Not a time-varying quantity
     708                self.quantities[q] = num.ravel(num.array(fin.variables[q], num.float)).reshape(M,3)
     709            else: # Time-varying, get the current timestep data
     710                self.quantities[q] = num.array(fin.variables[q][self.frame_number], num.float).reshape(M,3)
     711        fin.close()
     712        return self.quantities
     713
     714    def get_bounds(self):
     715        return [self.xmin, self.xmax, self.ymin, self.ymax]
     716
     717    def get_last_frame_number(self):
     718        return self.last_frame_number
     719
     720    def get_time(self):
     721        return self.time[self.frame_number]
    646722
    647723
     
    25302606            verbose = False,
    25312607            origin = None):
     2608   
    25322609    log.critical('sww2asc will soon be obsoleted - please use sww2dem')
    25332610    sww2dem(basename_in,
     
    25402617            verbose = verbose,
    25412618            origin = origin,
    2542         datum = 'WGS84',
    2543         format = 'asc')
     2619            datum = 'WGS84',
     2620            format = 'asc')
    25442621
    25452622
     
    63796456
    63806457
     6458
     6459
    63816460##
    63826461# @brief A class to write STS files.
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r7539 r7562  
    102102from anuga.geospatial_data.geospatial_data import ensure_geospatial
    103103
    104 from anuga.config import minimum_storable_height
    105 from anuga.config import minimum_allowed_height, maximum_allowed_speed
    106 from anuga.config import g, epsilon, beta_w, beta_w_dry,\
    107      beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
    108 from anuga.config import alpha_balance
    109 from anuga.config import optimise_dry_cells
    110 from anuga.config import optimised_gradient_limiter
    111 from anuga.config import use_edge_limiter
    112 from anuga.config import use_centroid_velocities
    113104from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    114105
     
    199190                                number_of_full_triangles=number_of_full_triangles)
    200191
     192        self.set_defaults()
     193
     194 
     195        self.forcing_terms.append(manning_friction_implicit)
     196        self.forcing_terms.append(gravity)
     197
     198        # Stored output
     199        self.store = True
     200        self.set_store_vertices_uniquely(False)
     201
     202        self.quantities_to_be_stored = {'elevation': 1,
     203                                        'stage': 2,
     204                                        'xmomentum': 2,
     205                                        'ymomentum': 2}
     206
     207
     208
     209
     210    ##
     211    # @brief Set default values, usually read in from a config file
     212    # @param flag
     213    def set_defaults(self):
     214        """Set the default values in this routine. That way we can inherit class
     215        and just over redefine the defaults for the new class
     216        """
     217
     218        from anuga.config import minimum_storable_height
     219        from anuga.config import minimum_allowed_height, maximum_allowed_speed
     220        from anuga.config import g, epsilon, beta_w, beta_w_dry,\
     221             beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters
     222        from anuga.config import alpha_balance
     223        from anuga.config import optimise_dry_cells
     224        from anuga.config import optimised_gradient_limiter
     225        from anuga.config import use_edge_limiter
     226        from anuga.config import use_centroid_velocities
     227
     228
     229
    201230        self.set_minimum_allowed_height(minimum_allowed_height)
    202 
    203231        self.maximum_allowed_speed = maximum_allowed_speed
     232
    204233        self.g = g
    205234        self.beta_w = beta_w
     
    214243        self.optimise_dry_cells = optimise_dry_cells
    215244
    216         self.use_new_mannings = False
    217         self.forcing_terms.append(manning_friction_implicit)
    218         self.forcing_terms.append(gravity)
    219 
    220         # Stored output
    221         self.store = True
    222         self.set_store_vertices_uniquely(False)
     245       
     246        self.set_new_mannings_function(False)
     247
    223248        self.minimum_storable_height = minimum_storable_height
    224         self.quantities_to_be_stored = {'elevation': 1,
    225                                         'stage': 2,
    226                                         'xmomentum': 2,
    227                                         'ymomentum': 2}
    228 
    229         # Limiters
     249
     250         # Limiters
    230251        self.use_edge_limiter = use_edge_limiter
    231252        self.optimised_gradient_limiter = optimised_gradient_limiter
    232         self.use_centroid_velocities = use_centroid_velocities
    233 
     253        self.use_centroid_velocities = use_centroid_velocities       
    234254
    235255    ##
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r7539 r7562  
    66776677                assert num.allclose(q, ref_flow)
    66786678
     6679        os.remove('Inflow_flowline_test.sww')
    66796680
    66806681    def test_volume_conservation_inflow(self):
     
    67766777
    67776778
     6779        os.remove('Inflow_volume_test.sww')
     6780
     6781
     6782       
    67786783    def test_volume_conservation_rain(self):
    67796784        """test_volume_conservation
     
    68816886            ref_volume += ys * delta_V
    68826887
     6888        os.remove('Rain_volume_test.sww')
    68836889
    68846890    def Xtest_rain_conservation_and_runoff(self):
     
    70847090        # Check that quantities have been stored correctly   
    70857091        from Scientific.IO.NetCDF import NetCDFFile
    7086         fid = NetCDFFile(domain.get_name() + '.sww')
     7092        sww_file = domain.get_name() + '.sww'
     7093        fid = NetCDFFile(sww_file)
    70877094
    70887095        x = fid.variables['x'][:]
     
    70917098        elevation = fid.variables['elevation'][:]
    70927099        fid.close()
     7100
     7101        os.remove(sww_file)
     7102       
    70937103                   
    70947104        assert len(stage.shape) == 2
     
    72617271                       % (q, ref_flow))
    72627272                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     7273
     7274        os.remove('inflow_flowline_test.sww')
    72637275
    72647276       
     
    75527564                               % (domain_depth, normal_depth))
    75537565
     7566        os.remove('Inflow_flowline_test.sww')
     7567
    75547568#################################################################################
    75557569
  • anuga_core/source/anuga/shallow_water_balanced/test_swb_distribute.py

    r7559 r7562  
    147147
    148148        # First order
    149         domain._order_ = 1
    150         domain.distribute_to_vertices_and_edges()
     149        domain.set_default_order(1)
     150        domain.distribute_to_vertices_and_edges()
     151       
    151152        assert num.allclose(L[1], val1)
    152153
    153154        # Second order
    154         domain._order_ = 2
    155         domain.beta_w = 0.9
    156         domain.beta_w_dry = 0.9
    157         domain.beta_uh = 0.9
    158         domain.beta_uh_dry = 0.9
    159         domain.beta_vh = 0.9
    160         domain.beta_vh_dry = 0.9
    161         domain.distribute_to_vertices_and_edges()
    162         assert num.allclose(L[1], [2.2, 4.9, 4.9])
     155        domain.set_default_order(2)
     156        domain.distribute_to_vertices_and_edges()
     157
     158        assert num.allclose(L[1], [1.00000000e-03,   5.99950000e+00,  5.99950000e+00])
     159
     160        assert num.allclose(val1, num.sum(L[1])/3)
    163161
    164162    def test_distribute_away_from_bed(self):
     
    192190        assert num.allclose(b[1], 0.0)
    193191
    194         domain._order_ = 1
    195         domain.distribute_to_vertices_and_edges()
     192        domain.set_default_order(1)
     193        domain.distribute_to_vertices_and_edges()
     194       
    196195        assert num.allclose(L[1], 1.77777778)
    197196
    198         domain._order_ = 2
    199         domain.beta_w = 0.9
    200         domain.beta_w_dry = 0.9
    201         domain.beta_uh = 0.9
    202         domain.beta_uh_dry = 0.9
    203         domain.beta_vh = 0.9
    204         domain.beta_vh_dry = 0.9
    205         domain.distribute_to_vertices_and_edges()
    206         assert num.allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
     197        domain.set_default_order(2)
     198        domain.distribute_to_vertices_and_edges()
     199
     200        assert num.allclose(L[1], [1.00000000e-03,   2.66616667e+00,   2.66616667e+00])
     201
     202        assert num.allclose(1.77777778, num.sum(L[1])/3)
    207203
    208204    def test_distribute_away_from_bed1(self):
     
    234230        assert num.allclose(b[1], 3.33333333)
    235231
    236         domain._order_ = 1
    237         domain.distribute_to_vertices_and_edges()
     232        domain.set_default_order(1)
     233        domain.distribute_to_vertices_and_edges()
     234       
    238235        assert num.allclose(L[1], 4.9382716)
    239236
    240         domain._order_ = 2
    241         domain.beta_w = 0.9
    242         domain.beta_w_dry = 0.9
    243         domain.beta_uh = 0.9
    244         domain.beta_uh_dry = 0.9
    245         domain.beta_vh = 0.9
    246         domain.beta_vh_dry = 0.9
    247         domain.distribute_to_vertices_and_edges()
    248         assert num.allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
     237        domain.set_default_order(2)
     238        domain.distribute_to_vertices_and_edges()
     239   
     240        assert num.allclose(L[1], [ 1.00000000e-03,   6.88207932e+00,   7.93173549e+00])
     241
    249242
    250243    def test_distribute_near_bed(self):
     
    283276                                domain.quantities['stage'].centroid_values[i])
    284277
    285         domain._order_ = 1
    286 
    287         domain.tight_slope_limiters = 0
    288         domain.distribute_to_vertices_and_edges()
    289         assert num.allclose(L[1], [0.1, 20.1, 20.1])
     278        domain.set_default_order(1)
     279        domain.distribute_to_vertices_and_edges()
     280
     281       
     282        assert num.allclose(L[1], [0.298,  20.001,  20.001])
    290283        for i in range(len(L)):
    291284            assert num.allclose(volumes[i], num.sum(L[i])/3)
    292285
    293         # Allow triangle to be flatter (closer to bed)
    294         domain.tight_slope_limiters = 1
    295 
    296         domain.distribute_to_vertices_and_edges()
    297         assert num.allclose(L[1], [0.298, 20.001, 20.001])
    298         for i in range(len(L)):
    299             assert num.allclose(volumes[i], num.sum(L[i])/3)
    300 
    301         domain._order_ = 2
    302 
    303         domain.tight_slope_limiters = 0
    304         domain.distribute_to_vertices_and_edges()
    305         assert num.allclose(L[1], [0.1, 20.1, 20.1])
    306         for i in range(len(L)):
    307             assert num.allclose(volumes[i], num.sum(L[i])/3)
    308 
    309         # Allow triangle to be flatter (closer to bed)
    310         domain.tight_slope_limiters = 1
    311 
    312         domain.distribute_to_vertices_and_edges()
    313         assert num.allclose(L[1], [0.298, 20.001, 20.001])
     286        domain.set_default_order(2)
     287        domain.distribute_to_vertices_and_edges()
     288
     289       
     290        assert num.allclose(L[1], [  0.1,  20.1,  20.1])
    314291        for i in range(len(L)):
    315292            assert num.allclose(volumes[i], num.sum(L[i])/3)
     
    378355
    379356        domain.distribute_to_vertices_and_edges()
    380         assert num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\
    381                num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\
    382                num.allclose(L[1], [4.19351461, 16.10548539, 20.001]) # old limiters
     357
     358
     359
     360
     361        assert num.allclose(L[1], [4.001, 16.15377472, 20.14522528],rtol=1.0e-3)
    383362       
    384363        for i in range(len(L)):
     
    440419        domain.distribute_to_vertices_and_edges()
    441420
    442         assert num.allclose(L[1,:], [-0.00825735775384,
    443                                      -0.00801881482869,
    444                                      0.0272445025825])
    445         assert num.allclose(X[1,:], [0.0143507718962,
    446                                      0.0142502147066,
    447                                      0.00931268339717])
    448         assert num.allclose(Y[1,:], [-0.000117062180693,
    449                                      7.94434448109e-005,
    450                                      -0.000151505429018])
     421
     422        assert num.allclose(L[1,:],[-0.01434766, -0.01292565, 0.03824164])
     423
     424        assert num.allclose(X[1,:],[ 0.01649334,  0.016267,    0.00515332])
     425
     426        assert num.allclose(Y[1,:],[-0.00027927,  0.00050728, -0.00041714])
    451427
    452428
  • anuga_core/source/anuga/shallow_water_balanced/test_swb_forcing_terms.py

    r7559 r7562  
    2424
    2525
    26 
    27 
    28 class Test_swb_clean(unittest.TestCase):
     26# Variable windfield implemented using functions
     27def speed(t, x, y):
     28    """Large speeds halfway between center and edges
     29
     30    Low speeds at center and edges
     31    """
     32
     33    from math import exp, cos, pi
     34
     35    x = num.array(x)
     36    y = num.array(y)
     37
     38    N = len(x)
     39    s = 0*x  #New array
     40
     41    for k in range(N):
     42        r = num.sqrt(x[k]**2 + y[k]**2)
     43        factor = exp(-(r-0.15)**2)
     44        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
     45
     46    return s
     47
     48def scalar_func(t, x, y):
     49    """Function that returns a scalar.
     50
     51    Used to test error message when numeric array is expected
     52    """
     53
     54    return 17.7
     55
     56def scalar_func_list(t, x, y):
     57    """Function that returns a scalar.
     58
     59    Used to test error message when numeric array is expected
     60    """
     61
     62    return [17.7]
     63
     64
     65def angle(t, x, y):
     66    """Rotating field
     67    """
     68    from math import atan, pi
     69
     70    x = num.array(x)
     71    y = num.array(y)
     72
     73    N = len(x)
     74    a = 0 * x    # New array
     75
     76    for k in range(N):
     77        r = num.sqrt(x[k]**2 + y[k]**2)
     78
     79        angle = atan(y[k]/x[k])
     80
     81        if x[k] < 0:
     82            angle += pi
     83
     84        # Take normal direction
     85        angle -= pi/2
     86
     87        # Ensure positive radians
     88        if angle < 0:
     89            angle += 2*pi
     90
     91        a[k] = angle/pi*180
     92
     93    return a
     94
     95
     96###############################################################################
     97
     98class Test_swb_forcing_terms(unittest.TestCase):
    2999    def setUp(self):
    30100        pass
     
    17231793
    17241794if __name__ == "__main__":
    1725     suite = unittest.makeSuite(Test_swb_clean, 'test')
     1795    suite = unittest.makeSuite(Test_swb_forcing_terms, 'test')
    17261796    runner = unittest.TextTestRunner(verbosity=1)
    17271797    runner.run(suite)
  • anuga_core/source/anuga/test_all.py

    r7318 r7562  
    2323
    2424# Directories that should not be searched for test files.
    25 exclude_dirs = ['pypar_dist',                        # Special requirements
    26                 '.svn',                              # subversion
     25exclude_dirs = ['pypar_dist',    # Special requirements
     26                '.svn',          # subversion
    2727                'props', 'wcprops', 'prop-base', 'text-base', 'tmp']
    2828
  • anuga_core/source/anuga_parallel/interface.py

    r7453 r7562  
    88from anuga_parallel.parallel_api import myid, numprocs, get_processor_name
    99from anuga_parallel.parallel_api import send, receive
    10 from anuga_parallel.parallel_api import pypar_available, barrier
     10from anuga_parallel.parallel_api import pypar_available, barrier, finalize
    1111
    1212from anuga_parallel.parallel_meshes import parallel_rectangle
  • anuga_core/source/anuga_parallel/test_parallel_distribute_domain.py

    r7505 r7562  
    3535
    3636
    37 from anuga_parallel.interface import distribute, myid, numprocs
     37from anuga_parallel.interface import distribute, myid, numprocs, finalize
    3838
    3939
     
    222222
    223223
     224
     225    finalize()
  • anuga_core/source/anuga_parallel/test_parallel_distribute_mesh.py

    r7507 r7562  
    1515from anuga_parallel.distribute_mesh import extract_hostmesh, rec_submesh, send_submesh
    1616
    17 from anuga_parallel.interface import myid, numprocs, barrier
     17from anuga_parallel.interface import myid, numprocs, barrier, finalize
    1818
    1919import numpy as num
     
    302302        distibute_three_processors()
    303303
    304 
     304    finalize()
    305305       
    306306
  • anuga_core/source/anuga_parallel/test_parallel_sw_flow.py

    r7505 r7562  
    2727
    2828
    29 from anuga_parallel.interface import distribute, myid, numprocs, send, receive, barrier
     29from anuga_parallel.interface import distribute, myid, numprocs, send, receive, barrier, finalize
    3030
    3131#--------------------------------------------------------------------------
     
    3737N = 25
    3838M = 25
    39 verbose = False
     39verbose = True
    4040
    4141#---------------------------------
     
    212212        evolution_test(parallel=True, G=G, seq_interpolation_points = interpolation_points)
    213213       
    214 
    215 
    216 
     214        finalize()
     215
     216
Note: See TracChangeset for help on using the changeset viewer.