Ignore:
Timestamp:
Feb 20, 2008, 4:01:00 PM (15 years ago)
Author:
sexton
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/MISG_2008/sensitivity_study.py

    r4937 r5066  
    33Water driven up a linear slope and time varying boundary,
    44similar to a beach environment
     5
     6Test case for MISG 2008, University of Wollongong
    57"""
    68
     
    1618from anuga.shallow_water import Time_boundary
    1719from anuga.shallow_water import Transmissive_boundary
     20from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
     21
     22from anuga.pmesh.mesh_interface import create_mesh_from_regions
    1823
    1924from Numeric import zeros, Float
     
    2530# Setup computational domain
    2631#------------------------------------------------------------------------------
    27 
    28 points, vertices, boundary = rectangular_cross(90, 30,
    29                                                len1=3.0,
    30                                                len2=1.0) # Basic mesh
    31 
    32 domain = Domain(points, vertices, boundary) # Create domain
     32y_extent=2000.0
     33bounding_polygon = [[0,0],[2000,0],[2000,2000],[0,2000]]
     34beach_polygon = [[5,900],[225,900],[225,1100],[5,1100]]
     35interior_regions = [[beach_polygon, 50]]
     36create_mesh_from_regions(bounding_polygon,
     37                         boundary_tags={'bottom':[0],
     38                                        'right':[1],
     39                                        'top':[2],
     40                                        'left':[3]},
     41                         maximum_triangle_area=2000,
     42                         interior_regions= interior_regions,
     43                         filename='misg.msh',
     44                         use_cache=False,
     45                         verbose=True)
     46domain = Domain('misg.msh',use_cache=False, verbose=True)
     47
     48
    3349domain.set_name('sensitivity')              # Output to file runup.sww
    3450domain.set_datadir('.')                     # Use current directory for output
    3551domain.tight_slope_limiters = 1
    3652domain.beta_h = 0
     53domain.set_minimum_storable_height(0.01)
     54
     55print 'number of elements: ', len(domain)
    3756
    3857#------------------------------------------------------------------------------
     
    4463    N = len(x)
    4564    z = zeros(N, Float)
    46 
    47 
     65   
    4866    for i in range(N):
    49         # linear bed slope           
    50         z[i] = -x[i]/15
    51 
    52         # IID noise
    53         z[i] += normal(0.0, 0.01)
     67        if x[i] <= (depth-c)/m:
     68            z[i] = m*x[i]+c
     69        else:
     70        # linear bed slope           
     71            #z[i] = slope*x[i] + depth-(slope)*(depth-c)/m
     72            z[i] = slope*x[i] + depth-(slope)*(depth-c)/m + 1.0*sin((x[i]-(depth-c)/m)/T) + 1.0*sin((y[i]-(depth-c)/m)/T)
     73
     74            # IID noise
     75            #z[i] += normal(-100.0, 5.0)
    5476       
    5577    return z
    56 
    57 domain.set_quantity('elevation', topography) # Use function for elevation
    58 domain.set_quantity('friction', 0.01)        # Constant friction
    59 domain.set_quantity('stage', -.4)            # Constant negative initial stage
    60 
    6178
    6279#------------------------------------------------------------------------------
     
    6986Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
    7087Bw = Time_boundary(domain=domain,     # Time dependent boundary 
    71                    f=lambda t: [(.1*sin(t*2*pi)-0.1) , 0.0, 0.0])
    72 
    73 # Associate boundary tags with boundary objects
    74 domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
    75 
    76 
    77 import time
    78 t0 = time.time()
    79 #------------------------------------------------------------------------------
    80 # Evolve system through time
    81 #------------------------------------------------------------------------------
    82 
    83 for t in domain.evolve(yieldstep = 0.1, finaltime = 50.0):
    84     domain.write_time()
    85 
    86 print 'That took %.2f seconds' %(time.time()-t0)
    87 
     88                   f=lambda t: [0.0, 0.0, 0.0])
     89
     90from Numeric import zeros, Float, arange
     91start=-1.0/40.0
     92finish=-1.8/40.0
     93n_vec=arange(start,finish,-2.0/400.0)
     94N=len(n_vec)
     95#N=1 # use for 2.5% and 5% change tests
     96max0=zeros(N, Float)
     97max1=zeros(N, Float)
     98max2=zeros(N, Float)
     99max3=zeros(N, Float)
     100mmax0=zeros(N, Float)
     101mmax1=zeros(N, Float)
     102mmax2=zeros(N, Float)
     103mmax3=zeros(N, Float)
     104output_file ='misg_influence_data.csv'
     105domain.store=True
     106fid_out = open(output_file, 'w')
     107s = 'variable, h0, h1, h2, h3, v0, v1, v2, v3\n'
     108fid_out.write(s)
     109
     110for i in range(N):
     111    manning = 0.01*1.0
     112    A=1.0*1.0
     113    w=1.0/100.0
     114    m=-1.0/20.0
     115    depth=-10.0
     116    c=10.0
     117    T=100.0*1.0
     118    #slope=-1.0/40.0*1.0
     119    slope = n_vec[i]
     120
     121    var = slope
     122    domain.set_quantity('elevation', topography) # Use function for elevation
     123    domain.set_quantity('friction', manning)
     124    domain.set_quantity('stage', 0.0)            # Constant negative initial stage
     125    domain.set_quantity('xmomentum', 0.0) # Use function for elevation
     126    domain.set_quantity('ymomentum', 0.0) # Use function for elevation
     127    swwname = 'sensitivity_2d_wavey_bottom'#+str(i)
     128    domain.set_name(swwname)
     129    Btd = Transmissive_Momentum_Set_Stage_boundary(domain=domain,
     130                                                   function=lambda t: [(0<t<pi/w)*A*sin(w*t), 0.0, 0.0])
     131    domain.set_boundary({'left': Br, 'right': Btd, 'top': Br, 'bottom': Br})
     132
     133    # locations to record wave height and velocity
     134    tri_id0 = domain.get_triangle_containing_point([-c/m,y_extent*0.5])
     135    tri_id1 = domain.get_triangle_containing_point([-c/m-10,y_extent*0.5])
     136    tri_id2 = domain.get_triangle_containing_point([-c/m-20,y_extent*0.5])
     137    tri_id3 = domain.get_triangle_containing_point([-c/m+50,y_extent*0.5])
     138     
     139    stage_centroid     = domain.quantities['stage'].centroid_values
     140    xmom_centroid      = domain.quantities['xmomentum'].centroid_values
     141    elevation_centroid = domain.quantities['elevation'].centroid_values
     142   
     143    elev0 = elevation_centroid[tri_id0]
     144    elev1 = elevation_centroid[tri_id1]
     145    elev2 = elevation_centroid[tri_id2]
     146    elev3 = elevation_centroid[tri_id3]
     147
     148    import time
     149    t0 = time.time()
     150    #------------------------------------------------------------------------------
     151    # Evolve system through time
     152    #------------------------------------------------------------------------------
     153    domain.time=0.0
     154    finaltime=4.0 #pi/w+150.0
     155    max=0
     156    for t in domain.evolve(yieldstep = 1.0, finaltime = finaltime):
     157
     158        if  stage_centroid[tri_id0]-elev0>max0[i]:
     159            max0[i]=stage_centroid[tri_id0]-elev0
     160        if  stage_centroid[tri_id1]-elev1>max1[i]:
     161            max1[i]=stage_centroid[tri_id1]-elev1
     162        if  stage_centroid[tri_id2]-elev2>max2[i]:
     163            max2[i]=stage_centroid[tri_id2]-elev2
     164        if  stage_centroid[tri_id3]-elev3>max3[i]:
     165            max3[i]=stage_centroid[tri_id3]-elev3
     166           
     167        if  abs(xmom_centroid[tri_id0])>mmax0[i]:
     168            mmax0[i]=abs(xmom_centroid[tri_id0])
     169        if  abs(xmom_centroid[tri_id1])>mmax1[i]:
     170            mmax1[i]=abs(xmom_centroid[tri_id1])
     171        if  abs(xmom_centroid[tri_id2])>mmax2[i]:
     172            mmax2[i]=abs(xmom_centroid[tri_id2])
     173        if  abs(xmom_centroid[tri_id3])>mmax3[i]:
     174            mmax3[i]=abs(xmom_centroid[tri_id3])
     175
     176    print 'finished cycle: ', i
     177    print 'That took %.2f seconds' %(time.time()-t0)
     178
     179
     180    s = '%.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f\n' %(var, max0[i], max1[i], max2[i], max3[i], mmax0[i], mmax1[i], mmax2[i], mmax3[i])
     181    fid_out.write(s)   
     182
     183from pylab import plot, ion, savefig, xlabel, ylabel, title, close, legend,hist,subplot
     184ion()
     185plot(n_vec,max2,'bo-',n_vec,max1,'go-',n_vec,max0,'ro-',n_vec,max3,'co-')
     186title('MISG testing: varying slope')
     187xlabel('slope')
     188ylabel('depth (m)')
     189#legend(['loc2','loc1','loc0','loc3'])
     190name='bed_slope'
     191savefig(name)
     192'''
     193#title('MISG testing: varying slope PDF')
     194subplot(2,2,1)
     195[h0,bins,patches]=hist(max0,10)
     196subplot(2,2,2)
     197[h1,bins,pa3ches]=hist(max1,10)
     198subplot(2,2,3)
     199[h2,bins,patches]=hist(max2,10)
     200subplot(2,2,4)
     201[h3,bins,patches]=hist(max3,10)
     202#plot(bins,h0,'bo-',bins,h1,'ro-',bins,h2,'go-',bins,h3,'co-')
     203xlabel('Runup')
     204ylabel('Frequency')
     205name='bed_slope_PDF'
     206savefig(name)
     207'''
     208
     209close('all')
Note: See TracChangeset for help on using the changeset viewer.