Changeset 3505


Ignore:
Timestamp:
Aug 17, 2006, 5:56:35 PM (18 years ago)
Author:
ole
Message:

Runup work

Files:
3 edited

Legend:

Unmodified
Added
Removed
  • documentation/user_manual/examples/runup.py

    r3486 r3505  
    4949#------------------------------------------------------------------------------
    5050
    51 from math import sin, pi
     51from math import sin, pi, exp
    5252Br = Reflective_boundary(domain)      # Solid reflective wall
    5353Bt = Transmissive_boundary(domain)    # Continue all values on boundary
    5454Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
    5555Bw = Time_boundary(domain=domain,     # Time dependent boundary 
    56                    f=lambda t: [(.1*sin(t*2*pi)-0.3), 0.0, 0.0])
     56                   f=lambda t: [(.1*sin(t*2*pi)-0.3) * exp(-t), 0.0, 0.0])
    5757
    5858# Associate boundary tags with boundary objects
  • documentation/user_manual/examples/runup_convergence.py

    r3487 r3505  
    22
    33Water driven up a linear slope and time varying boundary,
    4 similar to a beach environment
     4similar to a beach environment.
     5
     6The study area is discretised as a regular triangular grid 100m x 100m
    57"""
    68
     
    1012#------------------------------------------------------------------------------
    1113
     14import sys
     15
     16from pmesh.mesh_interface import create_mesh_from_regions
    1217from pyvolution.mesh_factory import rectangular_cross
    1318from pyvolution.shallow_water import Domain
     
    1520from pyvolution.shallow_water import Dirichlet_boundary
    1621from pyvolution.shallow_water import Time_boundary
    17 from pyvolution.shallow_water import Transmissive_boundary
     22from pyvolution.shallow_water import Transmissive_Momentum_Set_Stage_boundary
     23from pyvolution.util import file_function
     24from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure, axis
     25ion()
     26
     27
     28#------------------------------------------------------------------------------
     29# Model constants
     30
     31slope = -0.02     # 1:50 Slope
     32highest_point = 3 # Highest elevation (m)
     33sea_level = -2    # Mean sea level
     34wave_height = 3   # Solitary wave height H
     35simulation_name = 'runup_convergence'   
     36
     37
     38# Basin dimensions (m)
     39west = 0          # left boundary
     40east = 500        # right boundary
     41south = 0         # lower boundary
     42north = 100       # upper bourdary
    1843
    1944
     
    2247#------------------------------------------------------------------------------
    2348
    24 N = 10
    25 waveheight=5
    26 points, vertices, boundary = rectangular_cross(N, N, len1=100, len2=100)
     49# Structured mesh
     50dx = 10           # Resolution: Lenght of subdivisions on x axis (length)
     51dy = 10           # Resolution: Lenght of subdivisions on y axis (width)
     52
     53length = east-west
     54width = north-south
     55points, vertices, boundary = rectangular_cross(length/dx, width/dy, len1=length, len2=width,
     56                                               origin = (west, south))
    2757
    2858domain = Domain(points, vertices, boundary) # Create domain
    29 domain.set_name('runup_convergence')        # Output to bedslope.sww
     59
     60
     61
     62# Unstructured mesh
     63#polygon = [[east,north],[west,north],[west,south],[east,south]]
     64#meshname = simulation_name + '.msh'
     65#create_mesh_from_regions(polygon,
     66#                         boundary_tags={'top': [0], 'left': [1], 'bottom': [2], 'right': [3]},
     67#                         maximum_triangle_area=dx*dy/2,
     68#                         filename=meshname)
     69#                         #interior_regions=interior_regions)
     70#domain = Domain(meshname, use_cache=True, verbose = True)
     71
     72
     73
     74domain.set_name(simulation_name)
    3075
    3176
     
    3580
    3681def topography(x,y):
    37     slope = 0.02
    38     elevation  = 10
    39     return slope*x+elevation                      # linear bed slope
     82    return slope*x+highest_point  # Return linear bed slope bathymetry as vector
    4083
    4184domain.set_quantity('elevation', topography) # Use function for elevation
    42 domain.set_quantity('friction', 0. )         # Constant friction
    43 domain.set_quantity('stage', 0.0)            # Constant initial stage
     85domain.set_quantity('friction', 0.0 )        # Constant friction
     86domain.set_quantity('stage', sea_level)      # Constant initial stage
     87
    4488
    4589#------------------------------------------------------------------------------
     
    4791#------------------------------------------------------------------------------
    4892
    49 from math import sin, pi
     93from math import sin, pi, cosh, sqrt
    5094Br = Reflective_boundary(domain)      # Solid reflective wall
    51 Bt = Transmissive_boundary(domain)    # Continue all values on boundary
    5295Bd = Dirichlet_boundary([0.,0.,0.])   # Constant boundary values
     96
     97def waveform(t):
     98    return wave_height/cosh(t-5)**2
     99
    53100Bw = Time_boundary(domain=domain,     # Time dependent boundary 
    54                    f=lambda t: [((10<t<20)*waveheight), 0.0, 0.0])
     101                   f=lambda t: [waveform(t), 0.0, 0.0])
     102
     103# Time dependent boundary for stage, where momentum is set automatically
     104Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
     105
    55106
    56107# Associate boundary tags with boundary objects
    57 domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
     108domain.set_boundary({'left': Br, 'right': Bts, 'top': Br, 'bottom': Br})
    58109
    59110
     
    63114
    64115stagestep = []
    65 for t in domain.evolve(yieldstep = 10, finaltime = 100):
     116for t in domain.evolve(yieldstep = 1, finaltime = 100):
    66117    domain.write_time()
     118
    67119
    68120#-----------------------------------------------------------------------------
     
    70122#-----------------------------------------------------------------------------
    71123
     124
     125# Define line of gauges through center of domain
     126def gauge_line(west,east,north,south):
     127    from Numeric import arange
     128    x_vector = arange(west, (east-west)/2, 10) # Gauges every 1 meter from west to midpt
     129    y = (north+south)/2.
     130
     131    gauges = []
     132    for x in x_vector:
     133        gauges.append([x,y])
     134       
     135    return gauges, x_vector
     136
     137gauges, x_vector = gauge_line(west,east,north,south)
     138
     139# Obtain interpolated timeseries at gauges
     140f = file_function(domain.get_name()+'.sww',
     141                  quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
     142                  interpolation_points = gauges,
     143                  verbose = True,
     144                  use_cache = True)
     145
     146
     147# Find runup distance from western boundary through a linear search
     148max_stage = []
     149min_stage = []
     150runup_point = west
     151coastline = east       
     152for k, g in enumerate(gauges):
     153    z = f(0, point_id=k)[1] # Elevation
     154
     155    min_w = sys.maxint
     156    max_w = -min_w
     157    for i, t in enumerate(f.get_time()):
     158        w = f(t, point_id = k)[0]
     159        if w > max_w: max_w = w
     160        if w < min_w: min_w = w       
     161
     162    if max_w-z <= 0.01:  # Find first gauge where max depth > eps (runup)
     163        runup_point = g[0]
     164
     165    if min_w-z <= 0.01:  # Find first gauge where min depth > eps (coastline)
     166        coastline = g[0]       
     167       
     168    max_stage.append(max_w)
     169    min_stage.append(min_w)   
     170
     171
     172# Print
     173
     174print 'wave height [m]:                    ', wave_height
     175runup_height = topography(runup_point, (north+south)/2.)
     176print 'run up height [m]:                  ', runup_height
     177
     178runup_distance = runup_point-coastline
     179print 'run up distance from coastline [m]: ', runup_distance
     180
     181print 'Coastline (meters form west):       ', coastline
     182
     183
     184
     185# Store
     186filename = 'maxrunup'+str(wave_height)+'.csv'
     187fid = open(filename,'w')   
     188s = 'Depth at eastern boundary,Waveheight,Runup distance,Runup height\n'
     189fid.write(s)
     190
     191s = '%.2f,%.2f,%.2f\n' %(wave_height, runup_distance, runup_height)
     192fid.write(s)
     193
     194fid.close()
     195
     196# Plot
     197figure(1)
     198plot(x_vector, max_stage, 'g+',
     199     x_vector, min_stage, 'b+',     
     200     x_vector, topography(x_vector,(north+south)/2.), 'r-')
     201xlabel('x')
     202ylabel('stage')
     203title('Maximum stage for gauge line')
     204#axis([33000, 47000, -1000, 3000])
     205savefig('max_stage')
     206
     207close('all')
    72208   
    73209
  • inundation/pyvolution/shallow_water.py

    r3127 r3505  
    10591059        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
    10601060        value = self.function(self.domain.time)
    1061         q[0] = value[0]
     1061
     1062   
     1063        try:
     1064            # In case function returns a list of values
     1065            val = value[0]
     1066        except TypeError:
     1067            val = value
     1068
     1069        # Assign value to first component of q (stage)
     1070        try:
     1071            q[0] = float(val)
     1072        except:
     1073            msg = 'Supplied value could not be converted to float: val=', val     
     1074            raise Exception, msg   
     1075           
     1076           
     1077           
     1078       
    10621079        return q
    10631080
Note: See TracChangeset for help on using the changeset viewer.