| 1 | """Simple water flow example using ANUGA |
|---|
| 2 | |
|---|
| 3 | Water driven up a linear slope and time varying boundary, |
|---|
| 4 | similar to a beach environment |
|---|
| 5 | """ |
|---|
| 6 | |
|---|
| 7 | |
|---|
| 8 | #------------------------------------------------------------------------------ |
|---|
| 9 | # Import necessary modules |
|---|
| 10 | #------------------------------------------------------------------------------ |
|---|
| 11 | |
|---|
| 12 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
|---|
| 13 | from anuga.shallow_water import Domain |
|---|
| 14 | from anuga.shallow_water import Reflective_boundary |
|---|
| 15 | from anuga.shallow_water import Dirichlet_boundary |
|---|
| 16 | from anuga.shallow_water import Time_boundary |
|---|
| 17 | from anuga.shallow_water import Transmissive_boundary |
|---|
| 18 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
|---|
| 19 | from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files |
|---|
| 20 | from time import strftime, gmtime |
|---|
| 21 | from os import sep, environ, getenv, getcwd,umask |
|---|
| 22 | from anuga.utilities.polygon import Polygon_function |
|---|
| 23 | #------------------------------------------------------------------------------ |
|---|
| 24 | # Setup computational domain |
|---|
| 25 | #------------------------------------------------------------------------------ |
|---|
| 26 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
|---|
| 27 | |
|---|
| 28 | name ='steep_smflat_friction' |
|---|
| 29 | friction =[0, 0.02, 0.05, 0.1, 0.2, 0.3] |
|---|
| 30 | wave = [1, -1] #1 returns leading depression N-wave |
|---|
| 31 | #-1 returns leading crest N-wave |
|---|
| 32 | crest =[50, 450] |
|---|
| 33 | crestdepth = [2, -2, -6] |
|---|
| 34 | N = len (crest) |
|---|
| 35 | for i in range(N): |
|---|
| 36 | M = len (crestdepth) |
|---|
| 37 | for k in range (M): |
|---|
| 38 | O = len (friction) |
|---|
| 39 | for h in range (O): |
|---|
| 40 | B = len(wave) |
|---|
| 41 | for l in range(B): |
|---|
| 42 | length = (crest[i]+1000) |
|---|
| 43 | width = 20. |
|---|
| 44 | A = 1 |
|---|
| 45 | T = 2700 |
|---|
| 46 | umask(002) |
|---|
| 47 | time = strftime('%Y%m%d_%H%M%S',gmtime()) |
|---|
| 48 | |
|---|
| 49 | output_dir = sep+'d'+sep+'xrd'+sep+'gem'+sep+'5'+sep+'nhi'+sep+'inundation'+sep+'data'+sep+'idealised_bathymetry_study'+sep+'final_models'+sep+'friction'+sep+str(wave[l])+'_'+str(length)+'_'+str(A)+'_'+str(T)+'_'+str(crestdepth[k])+'_'+str(friction[h])+'_'+str(name)+sep |
|---|
| 50 | sww_file = str(name) |
|---|
| 51 | copy_code_files(output_dir,__file__,__file__) |
|---|
| 52 | |
|---|
| 53 | start_screen_catcher(output_dir) |
|---|
| 54 | dx = dy = .5 # Resolution: Length of subdivisions on both axes |
|---|
| 55 | boundary_polygon = [[0,0],[length,0],[length,width],[0,width]] |
|---|
| 56 | interior_polygon = [[140,0],[780+crest[i],0],[780+crest[i],20],[140,20]] |
|---|
| 57 | interior_polygon2 =[[790+crest[i],0],[999+crest[i],0],[999+crest[i],20],[790+crest[i],20]] |
|---|
| 58 | interior_regions = [[interior_polygon, 8], [interior_polygon2, 50]] |
|---|
| 59 | meshname = 'reef'+'.msh' |
|---|
| 60 | create_mesh_from_regions(boundary_polygon, |
|---|
| 61 | boundary_tags={'bottom': [0], |
|---|
| 62 | 'right': [1], |
|---|
| 63 | 'top': [2], |
|---|
| 64 | 'left': [3]}, |
|---|
| 65 | maximum_triangle_area=50, |
|---|
| 66 | filename=meshname, |
|---|
| 67 | interior_regions=interior_regions, |
|---|
| 68 | use_cache=False, |
|---|
| 69 | verbose=False) |
|---|
| 70 | |
|---|
| 71 | domain = Domain(meshname, use_cache=False, verbose=True) |
|---|
| 72 | |
|---|
| 73 | print 'Number of triangles = ', len(domain) |
|---|
| 74 | print 'The extent is ', domain.get_extent() |
|---|
| 75 | print domain.statistics() |
|---|
| 76 | |
|---|
| 77 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
|---|
| 78 | domain.set_minimum_storable_height(0.01) |
|---|
| 79 | domain.set_default_order(2) # Second order spatial approximation |
|---|
| 80 | domain.set_name(sww_file)# Output name |
|---|
| 81 | domain.set_datadir(output_dir) |
|---|
| 82 | |
|---|
| 83 | |
|---|
| 84 | #------------------------------------------------------------------------------ |
|---|
| 85 | # Setup initial conditions |
|---|
| 86 | #------------------------------------------------------------------------------ |
|---|
| 87 | |
|---|
| 88 | def topography(x,y): |
|---|
| 89 | """Complex topography defined by a function of vectors x and y |
|---|
| 90 | """ |
|---|
| 91 | |
|---|
| 92 | z =(-0.84*x)+(0.84*(720+crest[i]))+crestdepth[k]-4 |
|---|
| 93 | N = len (x) |
|---|
| 94 | for j in range(N): |
|---|
| 95 | |
|---|
| 96 | if x[j] < 180: |
|---|
| 97 | z[j] = -4+crestdepth[k] |
|---|
| 98 | |
|---|
| 99 | elif 179 < x[j] < 200: |
|---|
| 100 | z[j] = 0.2*x[j]-40+crestdepth[k] |
|---|
| 101 | #Reef Flat |
|---|
| 102 | elif 199 < x[j] < 700: |
|---|
| 103 | z[j] = -0.5+crestdepth[k] |
|---|
| 104 | ##Crest |
|---|
| 105 | elif 699 < x[j] < (700+crest[i]): |
|---|
| 106 | z[j] = +crestdepth[k] |
|---|
| 107 | |
|---|
| 108 | #Curve down |
|---|
| 109 | elif (699+crest[i]) < x[j] < (720+crest[i]): |
|---|
| 110 | z[j] = -0.01*(x[j]-(699+crest[i]))*(x[j]-(699+crest[i]))+crestdepth[k] |
|---|
| 111 | |
|---|
| 112 | return z |
|---|
| 113 | |
|---|
| 114 | |
|---|
| 115 | |
|---|
| 116 | domain.set_quantity('elevation', topography) # Use function for elevation |
|---|
| 117 | ## domain.set_quantity('friction', 0) # Constant friction |
|---|
| 118 | domain.set_quantity('friction', Polygon_function( [(boundary_polygon, 0.),(interior_polygon ,friction[h]), (interior_polygon2 , 0.)] ) )#changing friction over two polygons |
|---|
| 119 | domain.set_quantity('stage', 0.) # Constant negative initial stage |
|---|
| 120 | |
|---|
| 121 | |
|---|
| 122 | #------------------------------------------------------------------------------ |
|---|
| 123 | # Setup boundary conditions |
|---|
| 124 | #------------------------------------------------------------------------------ |
|---|
| 125 | |
|---|
| 126 | from math import sin, pi, exp, cos, cosh, sqrt |
|---|
| 127 | Br = Reflective_boundary(domain) # Solid reflective wall |
|---|
| 128 | Bt = Transmissive_boundary(domain) # Continue all values on boundary |
|---|
| 129 | Bd = Dirichlet_boundary([0.,0.,0.]) # Constant boundary values |
|---|
| 130 | Bw = Time_boundary(domain=domain, # Time dependent boundary |
|---|
| 131 | f=lambda t: [sin(2*pi*(t)/1010), -37, 0.0]) |
|---|
| 132 | g = 9.81 |
|---|
| 133 | offshore_depth = 250 |
|---|
| 134 | H_d_ratio = 0.008 |
|---|
| 135 | Xo = 67000 |
|---|
| 136 | po = 0.03 |
|---|
| 137 | def waveform(t): |
|---|
| 138 | return wave[l]*offshore_depth*(sqrt(g/offshore_depth)*t-Xo/offshore_depth)*sqrt(H_d_ratio*po)*H_d_ratio/cosh(sqrt(3*H_d_ratio*po/4)*(sqrt(g/offshore_depth)*t-Xo/offshore_depth))/cosh(sqrt(3*H_d_ratio*po/4)*(sqrt(g/offshore_depth)*t-Xo/offshore_depth)) |
|---|
| 139 | Bf = Transmissive_Momentum_Set_Stage_boundary(domain, waveform) |
|---|
| 140 | # Associate boundary tags with boundary objects |
|---|
| 141 | domain.set_boundary({'left': Bd, 'right': Bf, 'top': Br, 'bottom': Br}) |
|---|
| 142 | |
|---|
| 143 | |
|---|
| 144 | #------------------------------------------------------------------------------ |
|---|
| 145 | # Evolve system through time |
|---|
| 146 | #------------------------------------------------------------------------------ |
|---|
| 147 | |
|---|
| 148 | for t in domain.evolve(yieldstep = 5 , finaltime = 4500): |
|---|
| 149 | domain.write_time() |
|---|
| 150 | |
|---|
| 151 | """ |
|---|
| 152 | Generate time series of nominated "gauges" |
|---|
| 153 | Note, this script will only work if pylab is installed on the platform |
|---|
| 154 | |
|---|
| 155 | Inputs: |
|---|
| 156 | |
|---|
| 157 | production dirs: dictionary of production directories with a |
|---|
| 158 | association to that simulation run, eg high tide, |
|---|
| 159 | magnitude, etc. |
|---|
| 160 | |
|---|
| 161 | Outputs: |
|---|
| 162 | |
|---|
| 163 | * figures stored in same directory as sww file |
|---|
| 164 | * time series data stored in csv files in same directory as sww file |
|---|
| 165 | * elevation at nominated gauges (elev_output) |
|---|
| 166 | """ |
|---|
| 167 | |
|---|
| 168 | from os import getcwd, sep, altsep, mkdir, access, F_OK, remove |
|---|
| 169 | from anuga.abstract_2d_finite_volumes.util import sww2timeseries |
|---|
| 170 | |
|---|
| 171 | # nominate directory location of sww file with associated attribute |
|---|
| 172 | production_dirs = {output_dir: str(name)} |
|---|
| 173 | |
|---|
| 174 | # Generate figures |
|---|
| 175 | swwfiles = {} |
|---|
| 176 | for label_id in production_dirs.keys(): |
|---|
| 177 | file_loc = label_id |
|---|
| 178 | swwfile = file_loc + str(name)+'.sww' |
|---|
| 179 | swwfiles[swwfile] = label_id |
|---|
| 180 | print 'hello', swwfile |
|---|
| 181 | texname, elev_output = sww2timeseries(swwfiles, |
|---|
| 182 | sep+'d'+sep+'cit'+sep+'1'+sep+'cit'+sep+'natural_hazard_impacts'+sep+'inundation'+sep+'sandpits'+sep+'jbrowning'+sep+'anuga'+sep+'anuga_work'+sep+'development'+sep+'idealised_bathymetry_study'+sep+'final_models'+sep+'general'+sep+'gauges_smflat.csv', |
|---|
| 183 | production_dirs, |
|---|
| 184 | report = False, |
|---|
| 185 | reportname = '', |
|---|
| 186 | plot_quantity = ['stage', 'speed'], |
|---|
| 187 | generate_fig = False, |
|---|
| 188 | surface = False, |
|---|
| 189 | time_min = None, |
|---|
| 190 | time_max = None, |
|---|
| 191 | #time_unit = 'secs', |
|---|
| 192 | title_on = True, |
|---|
| 193 | verbose = True) |
|---|
| 194 | ## remove(output_dir+sep+str(name)+'.sww') |
|---|
| 195 | |
|---|
| 196 | |
|---|
| 197 | |
|---|
