Changeset 7869


Ignore:
Timestamp:
Jun 22, 2010, 5:52:28 PM (14 years ago)
Author:
hudson
Message:

Shallow water balanced tests fail, but no longer have errors.

Location:
trunk/anuga_core/source/anuga
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/__init__.py

    r7866 r7869  
    9393# Forcing
    9494#-----------------------------
    95 from anuga.shallow_water.forcing import Inflow, Rainfall
     95from anuga.shallow_water.forcing import Inflow, Rainfall, Wind_stress
    9696
    9797#-----------------------------

  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_basic.py

    r7866 r7869  
    133133
    134134    def test_rotate(self):
     135        from anuga.shallow_water.shallow_water_ext import rotate
    135136        normal = num.array([0.0, -1.0])
    136137
    … …  
    492493        # Setup boundary conditions
    493494        #--------------------------------------------------------------
    494         Br = Reflective_boundary(domain)                       # Reflective wall
    495         Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
     495        Br = anuga.Reflective_boundary(domain)                 # Reflective wall
     496        Bd = anuga.Dirichlet_boundary([final_runup_height, 0, 0])# Steady inflow
    496497
    497498        # All reflective to begin with (still water)
    … …  
    598599        from anuga.abstract_2d_finite_volumes.mesh_factory \
    599600                import rectangular_cross
    600         from data_manager import get_maximum_inundation_elevation
    601         from data_manager import get_maximum_inundation_location
    602         from data_manager import get_maximum_inundation_data
     601        from anuga.shallow_water.sww_interrogate import \
     602            get_maximum_inundation_elevation, get_maximum_inundation_location, \
     603            get_maximum_inundation_data
    603604
    604605        initial_runup_height = -0.4
    … …  
    632633        # Setup boundary conditions
    633634        #--------------------------------------------------------------
    634         Br = Reflective_boundary(domain)                       # Reflective wall
    635         Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
     635        Br = anuga.Reflective_boundary(domain)                 # Reflective wall
     636        Bd = anuga.Dirichlet_boundary([final_runup_height, 0, 0])# Steady inflow
    636637
    637638        # All reflective to begin with (still water)
    … …  
    787788        """
    788789
    789         from anuga.pmesh.mesh_interface import create_mesh_from_regions
    790         from anuga.abstract_2d_finite_volumes.mesh_factory \
    791                 import rectangular_cross
    792         from anuga.shallow_water import Domain
    793         from anuga.shallow_water import Reflective_boundary
    794         from anuga.shallow_water import Dirichlet_boundary
    795 
    796790        #-----------------------------------------------------------------
    797791        # Setup computational domain
    798792        #-----------------------------------------------------------------
    799         points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
    800         domain = Domain(points, vertices, boundary) # Create domain
     793        points, vertices, boundary = anuga.rectangular_cross(10, 10)# Basic mesh
     794        domain = anuga.Domain(points, vertices, boundary) # Create domain
    801795        domain.set_default_order(2)
    802796        domain.set_quantities_to_be_stored(None)
    … …  
    816810        # Setup boundary conditions
    817811        #----------------------------------------------------------------
    818         Br = Reflective_boundary(domain)           # Solid reflective wall
    819         Bd = Dirichlet_boundary([-0.2, 0., 0.])    # Constant boundary values
     812        Br = anuga.Reflective_boundary(domain)           # Solid reflective wall
     813        Bd = anuga.Dirichlet_boundary([-0.2, 0., 0.])    # Constant
    820814        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
    821815
    … …  
    879873        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    880874
    881         domain = Domain(points, vertices)
     875        domain = anuga.Domain(points, vertices)
    882876
    883877        #Set up for a gradient of (3,0) at mid triangle (bce)
    … …  
    898892        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
    899893
    900         domain.set_boundary({'exterior': Reflective_boundary(domain)})
     894        domain.set_boundary({'exterior': anuga.Reflective_boundary(domain)})
    901895
    902896        domain.optimise_dry_cells = True
    … …  
    916910
    917911    def test_second_order_flat_bed_onestep(self):
    918         from mesh_factory import rectangular
    919912
    920913        #Create basic mesh
    921         points, vertices, boundary = rectangular(6, 6)
     914        points, vertices, boundary = anuga.rectangular(6, 6)
    922915
    923916        #Create shallow water domain
    924         domain = Domain(points, vertices, boundary)
     917        domain = anuga.Domain(points, vertices, boundary)
    925918        domain.set_default_order(2)
    926919
    927920        # Boundary conditions
    928         Br = Reflective_boundary(domain)
    929         Bd = Dirichlet_boundary([0.1, 0., 0.])
     921        Br = anuga.Reflective_boundary(domain)
     922        Bd = anuga.Dirichlet_boundary([0.1, 0., 0.])
    930923        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
    931924
    … …  
    10271020
    10281021        # Create shallow water domain
    1029         domain = Domain(points, vertices, boundary)
     1022        domain = anuga.Domain(points, vertices, boundary)
    10301023        domain.smooth = False
    10311024        domain.default_order = 2
    10321025
    10331026        # Boundary conditions
    1034         Br = Reflective_boundary(domain)
    1035         Bd = Dirichlet_boundary([0.1, 0., 0.])
     1027        Br = anuga.Reflective_boundary(domain)
     1028        Bd = anuga.Dirichlet_boundary([0.1, 0., 0.])
    10361029        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
    10371030
    … …  
    10631056
    10641057        # Boundary conditions
    1065         Br = Reflective_boundary(domain)
    1066         Bd = Dirichlet_boundary([0.2, 0., 0.])
     1058        Br = anuga.Reflective_boundary(domain)
     1059        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
    10671060
    10681061        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
    … …  
    11011094
    11021095        # Boundary conditions
    1103         Br = Reflective_boundary(domain)
    1104         Bd = Dirichlet_boundary([0.2, 0., 0.])
     1096        Br = anuga.Reflective_boundary(domain)
     1097        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
    11051098
    11061099        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
    … …  
    11561149
    11571150        # Boundary conditions
    1158         Br = Reflective_boundary(domain)
    1159         Bd = Dirichlet_boundary([0.2, 0., 0.])
     1151        Br = anuga.Reflective_boundary(domain)
     1152        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
    11601153
    11611154        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
    … …  
    12011194
    12021195        # Boundary conditions
    1203         Br = Reflective_boundary(domain)
    1204         Bd = Dirichlet_boundary([0.2, 0., 0.])
     1196        Br = anuga.Reflective_boundary(domain)
     1197        Bd = anuga.Dirichlet_boundary([0.2, 0., 0.])
    12051198
    12061199        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
    … …  
    13481341
    13491342        # Boundary conditions
    1350         Br = Reflective_boundary(domain)
     1343        Br = anuga.Reflective_boundary(domain)
    13511344        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    13521345
    … …  
    15121505
    15131506        # Boundary conditions
    1514         Br = Reflective_boundary(domain)
     1507        Br = anuga.Reflective_boundary(domain)
    15151508        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    15161509
    … …  
    15561549        domain.set_quantity('elevation', Z)
    15571550
    1558         Br = Reflective_boundary(domain)
    1559         Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])
     1551        Br = anuga.Reflective_boundary(domain)
     1552        Bd = anuga.Dirichlet_boundary([inflow_stage, 0.0, 0.0])
    15601553        domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})
    15611554
    … …  
    16141607        import time, os
    16151608        from Scientific.IO.NetCDF import NetCDFFile
    1616         from data_manager import extent_sww
    16171609        from mesh_factory import rectangular_cross
    16181610
    … …  
    17281720
    17291721         tags = {}
    1730          b1 =  Dirichlet_boundary(conserved_quantities = num.array([0.0]))
    1731          b2 =  Dirichlet_boundary(conserved_quantities = num.array([1.0]))
    1732          b3 =  Dirichlet_boundary(conserved_quantities = num.array([2.0]))
     1722         b1 =  anuga.Dirichlet_boundary(conserved_quantities = num.array([0.0]))
     1723         b2 =  anuga.Dirichlet_boundary(conserved_quantities = num.array([1.0]))
     1724         b3 =  anuga.Dirichlet_boundary(conserved_quantities = num.array([2.0]))
    17331725         tags["1"] = b1
    17341726         tags["2"] = b2
    … …  
    19091901
    19101902        verbose = False
    1911 
    1912         from anuga.shallow_water import Domain
    1913         from anuga.pmesh.mesh_interface import create_mesh_from_regions
    1914         from anuga.geospatial_data.geospatial_data import Geospatial_data
    1915 
    19161903       
    19171904        # Get path where this test is run
    … …  
    19651952
    19661953        meshname = os.path.join(path, 'offending_mesh.msh')
    1967         create_mesh_from_regions(bounding_polygon,
     1954        anuga.create_mesh_from_regions(bounding_polygon,
    19681955                                 boundary_tags={'south': [0], 'east': [1],
    19691956                                                'north': [2], 'west': [3]},
    … …  
    19741961                                 verbose=verbose)
    19751962
    1976         domain = Domain(meshname, use_cache=False, verbose=verbose)
     1963        domain = anuga.Domain(meshname, use_cache=False, verbose=verbose)
    19771964
    19781965        #----------------------------------------------------------------------
    … …  
    20001987            os.remove(points_file)           
    20011988       
    2002         #finally:
    2003             # Cleanup regardless
    2004             #FIXME(Ole): Finally does not work like this in python2.4
    2005             #FIXME(Ole): Reinstate this when Python2.4 is out of the way
    2006             #FIXME(Ole): Python 2.6 apparently introduces something called 'with'           
    2007             #os.remove(meshname)
    2008             #os.remove(points_file)
    2009 
    20101989
    20111990    def test_fitting_example_that_crashed_2(self):
    … …  
    20222001
    20232002        verbose = False       
    2024 
    2025         from anuga.shallow_water import Domain
    2026         from anuga.pmesh.mesh_interface import create_mesh_from_regions
    2027         from anuga.geospatial_data import Geospatial_data
    20282003       
    20292004        # Get path where this test is run
    … …  
    20392014        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
    20402015
    2041         create_mesh_from_regions(bounding_polygon,
     2016        anuga.create_mesh_from_regions(bounding_polygon,
    20422017                                 boundary_tags={'south': [0],
    20432018                                                'east': [1],
    … …  
    20492024                                 verbose=verbose)
    20502025
    2051         domain = Domain(meshname, use_cache=True, verbose=verbose)
     2026        domain = anuga.Domain(meshname, use_cache=True, verbose=verbose)
    20522027       
    20532028        # Large test set revealed one problem
    … …  
    21312106        # Setup boundary conditions
    21322107        #------------------------------------------------------------------
    2133         Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
    2134         Br = Reflective_boundary(domain)              # Solid reflective wall
    2135         Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
     2108        Bi = anuga.Dirichlet_boundary([0.4, 0, 0])          # Inflow
     2109        Br = anuga.Reflective_boundary(domain)              # Solid reflective wall
     2110        Bo = anuga.Dirichlet_boundary([-5, 0, 0])           # Outflow
    21362111
    21372112        domain.set_boundary({'left': Bi, 'right': Bo, 'top': Br, 'bottom': Br})

  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py

    r7866 r7869  
    55from math import pi, sqrt
    66import tempfile
     7
     8import anuga
    79
    810from anuga.config import g, epsilon
    … …  
    6062                                          [30,30,30], [40,40,40]])
    6163
    62         D = Dirichlet_boundary([5, 2, 1])
    63         T = Transmissive_boundary(domain)
    64         R = Reflective_boundary(domain)
     64        D = anuga.Dirichlet_boundary([5, 2, 1])
     65        T = anuga.Transmissive_boundary(domain)
     66        R = anuga.Reflective_boundary(domain)
    6567        domain.set_boundary({'First': D, 'Second': T, 'Third': R})
    6668
    … …  
    140142                    (0, 1): 'Internal'}
    141143
    142         domain = Domain(points, vertices, boundary)
     144        domain = anuga.Domain(points, vertices, boundary)
    143145        domain.check_integrity()
    144146
    … …  
    150152                                          [30,30,30], [40,40,40]])
    151153
    152         D = Dirichlet_boundary([5, 2, 1])
    153         T = Transmissive_boundary(domain)
    154         R = Reflective_boundary(domain)
     154        D = anuga.Dirichlet_boundary([5, 2, 1])
     155        T = anuga.Transmissive_boundary(domain)
     156        R = anuga.Reflective_boundary(domain)
    155157        domain.set_boundary({'First': D, 'Second': T,
    156158                             'Third': R, 'Internal': None})
    … …  
    192194                                          [30,30,30], [40,40,40]])
    193195
    194         D = Dirichlet_boundary([5, 2, 1])
    195         T = Transmissive_stage_zero_momentum_boundary(domain)
    196         R = Reflective_boundary(domain)
     196        D = anuga.Dirichlet_boundary([5, 2, 1])
     197        T = anuga.Transmissive_stage_zero_momentum_boundary(domain)
     198        R = anuga.Reflective_boundary(domain)
    197199        domain.set_boundary({'First': D, 'Second': T, 'Third': R})
    198200
    … …  
    270272        #--------------------------------------------------------------
    271273        # Time dependent boundary
    272         Bt = Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0])
    273 
    274         Br = Reflective_boundary(domain)              # Reflective wall
     274        Bt = anuga.Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0])
     275
     276        Br = anuga.Reflective_boundary(domain)              # Reflective wall
    275277
    276278        domain.set_boundary({'left': Bt, 'right': Br, 'top': Br, 'bottom': Br})
    … …  
    326328
    327329        # Boundary conditions
    328         Br = Reflective_boundary(domain1)
    329         Bd = Dirichlet_boundary([0.3,0,0])
     330        Br = anuga.Reflective_boundary(domain1)
     331        Bd = anuga.Dirichlet_boundary([0.3,0,0])
    330332        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
    331333
    … …  
    370372
    371373        # Boundary conditions
    372         Br = Reflective_boundary(domain2)
    373         Bf = Field_boundary(domain1.get_name() + '.sww', domain2)
     374        Br = anuga.Reflective_boundary(domain2)
     375        Bf = anuga.Field_boundary(domain1.get_name() + '.sww', domain2)
    374376        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
    375377        domain2.check_integrity()
    … …  
    433435
    434436        # Boundary conditions
    435         Br = Reflective_boundary(domain1)
    436         Bd = Dirichlet_boundary([0.3,0,0])
     437        Br = anuga.Reflective_boundary(domain1)
     438        Bd = anuga.Dirichlet_boundary([0.3,0,0])
    437439        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
    438440
    … …  
    498500
    499501        # Boundary conditions
    500         Br = Reflective_boundary(domain2)
    501         Bf = Field_boundary(domain1.get_name() + '.sww',
     502        Br = anuga.Reflective_boundary(domain2)
     503        Bf = anuga.Field_boundary(domain1.get_name() + '.sww',
    502504                            domain2, verbose=False)
    503505        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
    … …  
    605607
    606608        # Boundary conditions
    607         Br = Reflective_boundary(domain1)
    608         Bd = Dirichlet_boundary([0.3, 0, 0])
     609        Br = anuga.Reflective_boundary(domain1)
     610        Bd = anuga.Dirichlet_boundary([0.3, 0, 0])
    609611        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
    610612
    … …  
    671673
    672674        # Boundary conditions
    673         Br = Reflective_boundary(domain2)
    674         Bf = Field_boundary(domain1.get_name() + '.sww',
     675        Br = anuga.Reflective_boundary(domain2)
     676        Bf = anuga.Field_boundary(domain1.get_name() + '.sww',
    675677                            domain2, mean_stage=mean_stage, verbose=False)
    676678
    … …  
    761763
    762764        # Create shallow water domain
    763         domain1 = Domain(points, vertices, boundary)
     765        domain1 = anuga.Domain(points, vertices, boundary)
    764766
    765767        domain1.reduction = mean
    … …  
    784786
    785787        # Boundary conditions
    786         Br = Reflective_boundary(domain1)
    787         Bd = Dirichlet_boundary([0.3, 0, 0])
     788        Br = anuga.Reflective_boundary(domain1)
     789        Bd = anuga.Dirichlet_boundary([0.3, 0, 0])
    788790        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
    789791
    … …  
    846848
    847849        # Boundary conditions
    848         Br = Reflective_boundary(domain2)
    849         Bf = Field_boundary(domain1.get_name() + '.sww',
     850        Br = anuga.Reflective_boundary(domain2)
     851        Bf = anuga.Field_boundary(domain1.get_name() + '.sww',
    850852                            domain2, mean_stage=1, verbose=False)
    851853

  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_forcing_terms.py

    r7866 r7869  
    552552        os.remove(filename + '.tms')
    553553
    554         W = Wind_stress(F)
     554        W = anuga.Wind_stress(F)
    555555
    556556        domain.forcing_terms = []
    … …  
    11631163        # on a circle affecting triangles #0 and #1 (bac and bce)
    11641164        domain.forcing_terms = []
    1165         I = Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1)
     1165        I = anuga.Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1)
    11661166        domain.forcing_terms.append(I)
    11671167       
    … …  
    12901290        # for a range of stage values
    12911291        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
    1292             print stage
    1293 
    12941292            domain.time = 0.0
    12951293            domain.set_quantity('stage', stage)
    … …  
    13021300            for t in domain.evolve(yieldstep=dt, finaltime=5.0):
    13031301                volume = domain.quantities['stage'].get_integral()
    1304                 print t, volume, predicted_volume
    13051302                assert num.allclose (volume, predicted_volume)
    13061303                predicted_volume = predicted_volume - 2.0/pi/100/dt # Why 100?
    … …  
    13091306        # range of stage values
    13101307        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
    1311             print stage
    1312 
    13131308            domain.time = 0.0
    13141309            domain.set_quantity('stage', stage)
    … …  
    13231318            for t in domain.evolve(yieldstep=dt, finaltime=5.0):
    13241319                volume = domain.quantities['stage'].get_integral()
    1325 
    1326                 print t, volume
    13271320                assert num.allclose(volume, initial_volume)
    13281321
    … …  
    14031396                    domain.forcing_terms.append(fixed_inflow)
    14041397               
    1405                 ref_flow = fixed_inflow.rate*number_of_inflowsg
     1398                ref_flow = fixed_inflow.rate*number_of_inflows
    14061399
    14071400                # Compute normal depth on plane using Mannings equation
    … …  
    14171410                #--------------------------------------------------------------
    14181411
    1419                 Br = Reflective_boundary(domain)
     1412                Br = anuga.Reflective_boundary(domain)
    14201413               
    14211414                # Define downstream boundary based on predicted depth
    … …  
    14231416                    return (-slope*length) + normal_depth
    14241417               
    1425                 Bt = Transmissive_momentum_set_stage_boundary(domain=domain,
    1426                                                               function=normal_depth_stage_downstream)
    1427                
    1428 
    1429                
     1418                Bt = anuga.Transmissive_momentum_set_stage_boundary(
     1419                        domain=domain, function=normal_depth_stage_downstream)
    14301420
    14311421                domain.set_boundary({'left': Br,
    … …  
    14331423                                     'top': Br,
    14341424                                     'bottom': Br})
    1435 
    14361425
    14371426
    … …  
    16171606
    16181607       
    1619        
    1620     def Xtest_friction_dependent_flow_using_flowline(self):
    1621         """test_friction_dependent_flow_using_flowline
    1622        
    1623         Test the internal flow (using flowline) as a function of
    1624         different values of Mannings n and different slopes.
    1625        
    1626         Flow is applied in the form of boundary conditions with fixed momentum.
    1627         """
    1628 
    1629         verbose = True
    1630 
    1631         #----------------------------------------------------------------------
    1632         # Import necessary modules
    1633         #----------------------------------------------------------------------
    1634 
    1635         from anuga.abstract_2d_finite_volumes.mesh_factory \
    1636                 import rectangular_cross
    1637         from anuga.shallow_water import Domain
    1638         from anuga.shallow_water.shallow_water_domain import Reflective_boundary
    1639         from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
    1640         from anuga.shallow_water.forcing import Inflow
    1641         from anuga.shallow_water.data_manager \
    1642                 import get_flow_through_cross_section
    1643         from anuga.abstract_2d_finite_volumes.util \
    1644                 import sww2csv_gauges, csv2timeseries_graphs
    1645 
    1646 
    1647         #----------------------------------------------------------------------
    1648         # Setup computational domain
    1649         #----------------------------------------------------------------------
    1650 
    1651         finaltime = 1000.0
    1652 
    1653         length = 300.
    1654         width  = 20.
    1655         dx = dy = 5       # Resolution: of grid on both axes
    1656        
    1657         # Input parameters
    1658         uh = 1.0
    1659         vh = 0.0
    1660         d = 1.0
    1661        
    1662         ref_flow = uh*d*width # 20 m^3/s in the x direction across entire domain
    1663 
    1664         points, vertices, boundary = rectangular_cross(int(length/dx),
    1665                                                        int(width/dy),
    1666                                                        len1=length,
    1667                                                        len2=width)
    1668 
    1669         for mannings_n in [0.035]:          #[0.0, 0.012, 0.035]:
    1670             for slope in [1.0/300]:         #[0.0, 1.0/300, 1.0/150]:
    1671                 # Loop over a range of bedslopes representing
    1672                 # sub to super critical flows
    1673                 if verbose:
    1674                     print
    1675                     print 'Slope:', slope, 'Mannings n:', mannings_n
    1676                 domain = Domain(points, vertices, boundary)   
    1677                 domain.set_name('Inflow_flowline_test')     # Output name
    1678 
    1679                 #--------------------------------------------------------------
    1680                 # Setup initial conditions
    1681                 #--------------------------------------------------------------
    1682 
    1683                 def topography(x, y):
    1684                     z = -x * slope
    1685                     return z
    1686 
    1687                 # Use function for elevation
    1688                 domain.set_quantity('elevation', topography)
    1689                 # Constant friction
    1690                 domain.set_quantity('friction', mannings_n)
    1691                
    1692                 #domain.set_quantity('stage', expression='elevation')
    1693                      
    1694                 # Set initial flow as depth=1m, uh=1.0 m/s, vh = 0.0
    1695                 # making it 20 m^3/s across entire domain
    1696                 domain.set_quantity('stage', expression='elevation + %f' % d)
    1697                 domain.set_quantity('xmomentum', uh)
    1698                 domain.set_quantity('ymomentum', vh)               
    1699 
    1700                 #--------------------------------------------------------------
    1701                 # Setup boundary conditions
    1702                 #--------------------------------------------------------------
    1703 
    1704                 Br = Reflective_boundary(domain)      # Solid reflective wall
    1705                
    1706                 # Constant flow in and out of domain
    1707                 # Depth = 1m, uh=1 m/s, i.e. a flow of 20 m^3/s
    1708                 # across boundaries
    1709                 Bi = Dirichlet_boundary([d, uh, vh])
    1710                 Bo = Dirichlet_boundary([-length*slope+d, uh, vh])
    1711                 #Bo = Dirichlet_boundary([-100, 0, 0])
    1712 
    1713                 domain.set_boundary({'left': Bi, 'right': Bo,
    1714                                      'top': Br,  'bottom': Br})
    1715 
    1716                 #--------------------------------------------------------------
    1717                 # Evolve system through time
    1718                 #--------------------------------------------------------------
    1719 
    1720                 for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
    1721                     if verbose :
    1722                         print domain.timestepping_statistics()
    1723                         print domain.volumetric_balance_statistics()
    1724 
    1725                 # 90 degree flowline at 200m
    1726                 q = domain.get_flow_through_cross_section([[200.0,  0.0],
    1727                                                            [200.0, 20.0]])
    1728                 msg = ('Predicted flow was %f, should have been %f'
    1729                        % (q, ref_flow))
    1730                 if verbose:
    1731                     print ('90 degree flowline: ANUGA = %f, Ref = %f'
    1732                            % (q, ref_flow))
    1733 
    1734                 # 45 degree flowline at 200m
    1735                 q = domain.get_flow_through_cross_section([[200.0,  0.0],
    1736                                                            [220.0, 20.0]])
    1737                 msg = ('Predicted flow was %f, should have been %f'
    1738                        % (q, ref_flow))
    1739                 if verbose:
    1740                     print ('45 degree flowline: ANUGA = %f, Ref = %f'
    1741                            % (q, ref_flow))
    1742 
    1743                 # Stage recorder (gauge) in middle of plane at 200m
    1744                 x = 200.0
    1745                 y = 10.00               
    1746                 w = domain.get_quantity('stage').\
    1747                         get_values(interpolation_points=[[x, y]])[0]
    1748                 z = domain.get_quantity('elevation').\
    1749                         get_values(interpolation_points=[[x, y]])[0]
    1750                 domain_depth = w-z
    1751 
    1752                 xmom = domain.get_quantity('xmomentum').\
    1753                         get_values(interpolation_points=[[x, y]])[0]
    1754                 ymom = domain.get_quantity('ymomentum').\
    1755                         get_values(interpolation_points=[[x, y]])[0]           
    1756                 if verbose:                   
    1757                     print ('At interpolation point (h, uh, vh): ',
    1758                            domain_depth, xmom, ymom)
    1759                     print 'uh * d * width = ', xmom*domain_depth*width
    1760                                
    1761                 if slope > 0.0:
    1762                     # Compute normal depth at gauge location using Manning eqn
    1763                     # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
    1764                     normal_depth = (ref_flow*mannings_n/(slope**0.5*width))**0.6
    1765                     if verbose:
    1766                         print ('Depth: ANUGA = %f, Mannings = %f'
    1767                                % (domain_depth, normal_depth))
    1768 
    1769         os.remove('Inflow_flowline_test.sww')
    1770 
    1771 
    1772 #################################################################################
     1608 ################################################################################
    17731609
    17741610if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.