source: trunk/anuga_work/development/gareth/tests/runup/runup.py @ 8308

Last change on this file since 8308 was 8308, checked in by davies, 12 years ago

Added extrapolate_velocity_second_order option in config.py
and getting various tests working in gareth/tests

File size: 4.2 KB
Line 
1"""Runup example from the manual, slightly modified
2"""
3#---------
4#Import Modules
5#--------
6from sys import path
7import anuga
8
9import numpy
10
11import struct
12#import scipy
13
14#from Numeric import *
15
16from math import sin, pi, exp
17#from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
18#from anuga.shallow_water.shallow_water_domain import Domain as Domain
19#from shallow_water_balanced_steve.swb_domain import *
20#import shallow_water_balanced_steve.swb_domain
21#from shallow_water_balanced_steve.swb_domain import Domain as Domain
22#path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/shallow_water_balanced_steve')
23#from swb_domain import *
24#path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
25from balanced_basic import *
26#---------
27#Setup computational domain
28#---------
29points, vertices, boundary = anuga.rectangular_cross(40,40)
30
31domain=Domain(points,vertices,boundary)    # Create Domain
32domain.set_name('runup_v2')                         # Output to file runup.sww
33domain.set_datadir('.')                          # Use current folder
34domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
35domain.set_store_vertices_uniquely(True)
36#------------------
37# Define topography
38#------------------
39
40def topography(x,y):
41        return -x/2 #+0.05*numpy.sin((x+y)*50.0) #+0.1*(numpy.random.rand(len(x)) -0.5)       # Linear bed slope + small random perturbation
42
43def stagefun(x,y):
44    #stg=-0.2*(x<0.5) -0.1*(x>=0.5)
45    stg=-0.2 # Stage
46    #topo=topography(x,y) #Bed
47    return stg #*(stg>topo) + topo*(stg<=topo)
48
49domain.set_quantity('elevation',topography)     # Use function for elevation
50domain.get_quantity('elevation').smooth_vertex_values() # Steve's fix -- without this, substantial artificial velcities are generated everywhere in the domain. With this fix, there are artificial velocities near the coast, but not elsewhere.
51
52domain.set_quantity('friction',0.00)             # Constant friction
53
54domain.set_quantity('stage', stagefun)              # Constant negative initial stage
55
56#--------------------------
57# Setup boundary conditions
58#--------------------------
59Br=anuga.Reflective_boundary(domain)            # Solid reflective wall
60Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
61Bd=anuga.Dirichlet_boundary([-0.2,0.,0.])       # Constant boundary values -- not used in this example
62Bw=anuga.Time_boundary(domain=domain,
63        f=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,0.0,0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.
64
65#----------------------------------------------
66# Associate boundary tags with boundary objects
67#----------------------------------------------
68domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
69
70#------------------------------
71#Evolve the system through time
72#------------------------------
73
74#xwrite=open("xvel.out","wb")
75#ywrite=open("yvel.out","wb")
76## Set print options to be compatible with file writing via the 'print' statement
77#numpy.set_printoptions(threshold=numpy.nan, linewidth=numpy.nan)
78
79for t in domain.evolve(yieldstep=0.2,finaltime=30.0):
80    print domain.timestepping_statistics()
81
82#    momx=domain.quantities['xmomentum'].centroid_values
83#    momy=domain.quantities['ymomentum'].centroid_values
84#    #mom2=domain.quantities['xmomentum'].vertex_values
85#    dep1=(domain.quantities['stage'].centroid_values-domain.quantities['elevation'].centroid_values+1.0e-06)
86#    #dep2=(domain.quantities['stage'].vertex_values-domain.quantities['elevation'].vertex_values+1.0e-06)
87#    velx=momx/dep1
88#    vely=momy/dep1
89#    #vel2=mom2/dep2
90#    #print vel1.max(), vel2.max()
91#    #print vel1.min(), vel2.min()
92#    #xwrite.write(velx)
93#    #print >> xwrite, str(velx)
94#    #print >> ywrite, str(vely)
95#    for i in velx:
96#        data=struct.pack('f',i)
97#        xwrite.write(data)
98#
99#    for j in vely:
100#        data2=struct.pack('f',j)
101#        ywrite.write(data2)
102#
103#    #print >> xwrite, \n
104#    #numpy.savetxt("xvel.txt",velx)
105#    #numpy.savetxt("yvel.txt",vely)
106#    #velx.tofile(xwrite," ")
107#    #vely.tofile(ywrite," ")
108#
109#xwrite.close()
110#ywrite.close()
111
112print 'Finished'
Note: See TracBrowser for help on using the repository browser.