[5093] | 1 | """ |
---|
| 2 | |
---|
| 3 | Script for running a breaking wave simulation of Jon Hinwoods wave tank. |
---|
| 4 | Note: this is based on the frinction_ua_flume_2006 structure. |
---|
| 5 | |
---|
| 6 | |
---|
| 7 | Duncan Gray, GA - 2007 |
---|
| 8 | |
---|
| 9 | |
---|
| 10 | |
---|
| 11 | """ |
---|
| 12 | |
---|
| 13 | |
---|
| 14 | #---------------------------------------------------------------------------- |
---|
| 15 | # Import necessary modules |
---|
| 16 | #---------------------------------------------------------------------------- |
---|
| 17 | |
---|
| 18 | # Standard modules |
---|
| 19 | import time |
---|
| 20 | from time import localtime, strftime |
---|
| 21 | import sys |
---|
| 22 | from shutil import copy |
---|
| 23 | from os import path, sep |
---|
| 24 | from os.path import dirname #, basename |
---|
[5095] | 25 | from math import sin, pi |
---|
[5093] | 26 | |
---|
| 27 | # Related major packages |
---|
| 28 | from anuga.shallow_water import Domain, Reflective_boundary, \ |
---|
| 29 | Dirichlet_boundary, Time_boundary, \ |
---|
| 30 | File_boundary, \ |
---|
| 31 | Transmissive_Momentum_Set_Stage_boundary |
---|
| 32 | from anuga.fit_interpolate.interpolate import interpolate_sww2csv |
---|
[5095] | 33 | from anuga.abstract_2d_finite_volumes.util import copy_code_files, \ |
---|
| 34 | file_function |
---|
[5190] | 35 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
| 36 | |
---|
[5095] | 37 | |
---|
| 38 | from anuga.shallow_water.data_manager import start_screen_catcher |
---|
| 39 | |
---|
[5093] | 40 | from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\ |
---|
| 41 | import File_boundary_time |
---|
| 42 | |
---|
| 43 | # Scenario specific imports |
---|
| 44 | import project # Definition of file names and polygons |
---|
| 45 | import create_mesh |
---|
| 46 | |
---|
| 47 | def elevation_function(x,y): |
---|
| 48 | from Numeric import zeros, size, Float |
---|
| 49 | |
---|
| 50 | z = zeros(size(x), Float) |
---|
| 51 | for i in range(len(x)): |
---|
| 52 | if x[i] <= 0: |
---|
| 53 | # swash |
---|
| 54 | z[i] = -0.05*x[i] |
---|
[5154] | 55 | #z[i] = -30 |
---|
[5093] | 56 | elif x[i] <= 1000: |
---|
| 57 | # surf |
---|
| 58 | z[i] = -0.15*(x[i]**(0.5)) |
---|
[5154] | 59 | #z[i] = -30 |
---|
[5093] | 60 | else: |
---|
| 61 | # Nominal shoreface |
---|
| 62 | z[i] = -0.01888*(x[i]**(0.8)) |
---|
[5154] | 63 | #z[i] = -30 |
---|
[5093] | 64 | #>>> -0.15*(1000**0.5) |
---|
| 65 | #-4.7434164902525691 |
---|
| 66 | |
---|
| 67 | #>>> 0.01888*(1000**0.8) |
---|
| 68 | #4.7424415826900885 |
---|
[5095] | 69 | |
---|
| 70 | #>>> -0.01888*10000**0.8 |
---|
| 71 | #-29.922783473665834 |
---|
[5093] | 72 | return z |
---|
| 73 | |
---|
[5190] | 74 | def main(friction=0.01, outputdir_name=None, |
---|
| 75 | is_structured=False, |
---|
| 76 | is_trial_run=False): |
---|
| 77 | |
---|
| 78 | # Structured mesh variables |
---|
[5192] | 79 | dx = 1.0 |
---|
[5190] | 80 | dy = dx |
---|
| 81 | L = 2080. |
---|
| 82 | W = dx |
---|
| 83 | outputdir_name += '_dx' + str(dx) |
---|
| 84 | |
---|
| 85 | basename = 'dx_' + str(dx) |
---|
[5093] | 86 | if is_trial_run is True: |
---|
| 87 | outputdir_name += '_test' |
---|
[5154] | 88 | yieldstep = 1 |
---|
[5192] | 89 | finaltime = 200. |
---|
[5095] | 90 | maximum_triangle_area=3000 |
---|
| 91 | thinner=True |
---|
[5093] | 92 | else: |
---|
[5154] | 93 | yieldstep = 1. |
---|
[5192] | 94 | finaltime = 1700 #1600 |
---|
[5154] | 95 | maximum_triangle_area = 100 |
---|
[5095] | 96 | thinner=True |
---|
[5093] | 97 | |
---|
| 98 | |
---|
| 99 | pro_instance = project.Project(['data','near_shore_PMD'], |
---|
[5095] | 100 | outputdir_name=outputdir_name) |
---|
[5093] | 101 | print "The output dir is", pro_instance.outputdir |
---|
| 102 | copy_code_files(pro_instance.outputdir,__file__, |
---|
| 103 | dirname(project.__file__) \ |
---|
| 104 | + sep + project.__name__+'.py') |
---|
[5095] | 105 | copy (pro_instance.codedir + 'run_beach.py', |
---|
| 106 | pro_instance.outputdir + 'run_beach.py') |
---|
[5093] | 107 | copy (pro_instance.codedir + 'create_mesh.py', |
---|
| 108 | pro_instance.outputdir + 'create_mesh.py') |
---|
| 109 | |
---|
[5154] | 110 | #mesh_filename = pro_instance.meshdir + basename + '.tsh' |
---|
[5093] | 111 | mesh_filename = pro_instance.meshdir + basename + '.msh' |
---|
| 112 | |
---|
| 113 | #-------------------------------------------------------------------------- |
---|
| 114 | # Copy scripts to output directory and capture screen |
---|
| 115 | # output to file |
---|
| 116 | #-------------------------------------------------------------------------- |
---|
| 117 | |
---|
| 118 | # creates copy of code in output dir |
---|
| 119 | if is_trial_run is False: |
---|
[5095] | 120 | start_screen_catcher(pro_instance.outputdir) #, rank, pypar.size()) |
---|
[5093] | 121 | |
---|
| 122 | print 'USER: ', pro_instance.user |
---|
| 123 | #------------------------------------------------------------------------- |
---|
| 124 | # Create the triangular mesh |
---|
| 125 | #------------------------------------------------------------------------- |
---|
| 126 | |
---|
| 127 | # this creates the mesh |
---|
| 128 | #gate_position = 12.0 |
---|
[5190] | 129 | if is_structured is True: |
---|
[5093] | 130 | |
---|
[5190] | 131 | # Original struc mesh values |
---|
| 132 | # dx = 10. |
---|
| 133 | # dy = dx |
---|
| 134 | # L = 2080. |
---|
| 135 | # W = dx |
---|
| 136 | |
---|
| 137 | # structured mesh |
---|
| 138 | points, vertices, boundary = rectangular_cross( |
---|
| 139 | int(L/dx), int(W/dy), L, W, (-80.0, -W*0.5)) |
---|
| 140 | |
---|
| 141 | domain = Domain(points, vertices, boundary) |
---|
| 142 | else: |
---|
| 143 | create_mesh.generate(mesh_filename, |
---|
| 144 | maximum_triangle_area=maximum_triangle_area, |
---|
| 145 | thinner=thinner) |
---|
| 146 | |
---|
| 147 | head,tail = path.split(mesh_filename) |
---|
| 148 | copy (mesh_filename, |
---|
| 149 | pro_instance.outputdir + tail ) |
---|
| 150 | domain = Domain(mesh_filename, use_cache = False, verbose = True) |
---|
| 151 | domain = Domain(points, vertices, boundary) |
---|
[5093] | 152 | #------------------------------------------------------------------------- |
---|
| 153 | # Setup computational domain |
---|
| 154 | #------------------------------------------------------------------------- |
---|
| 155 | |
---|
| 156 | |
---|
| 157 | print 'Number of triangles = ', len(domain) |
---|
| 158 | print 'The extent is ', domain.get_extent() |
---|
| 159 | print domain.statistics() |
---|
| 160 | |
---|
| 161 | |
---|
| 162 | domain.set_name(basename) |
---|
| 163 | domain.set_datadir(pro_instance.outputdir) |
---|
[5154] | 164 | domain.set_default_order(2) # Use second order spatial scheme |
---|
[5190] | 165 | domain.set_timestepping_method('rk2') |
---|
| 166 | domain.set_default_order(2) |
---|
[5191] | 167 | domain.use_edge_limiter = True |
---|
| 168 | domain.tight_slope_limiters = True |
---|
[5190] | 169 | |
---|
[5191] | 170 | domain.beta_w = 0.6 |
---|
| 171 | domain.beta_uh = 0.6 |
---|
| 172 | domain.beta_vh = 0.6 |
---|
| 173 | domain.beta_h = 0.0 |
---|
[5190] | 174 | |
---|
[5093] | 175 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
| 176 | domain.set_minimum_storable_height(0.001) |
---|
| 177 | #domain.set_store_vertices_uniquely(True) # for writting to sww |
---|
| 178 | |
---|
| 179 | #------------------------------------------------------------------------- |
---|
| 180 | # Setup initial conditions |
---|
| 181 | #------------------------------------------------------------------------- |
---|
| 182 | |
---|
| 183 | domain.set_quantity('stage', 0.0) |
---|
| 184 | domain.set_quantity('friction', friction) |
---|
| 185 | domain.set_quantity('elevation', elevation_function) |
---|
| 186 | |
---|
| 187 | |
---|
| 188 | print 'Available boundary tags', domain.get_boundary_tags() |
---|
| 189 | |
---|
| 190 | # Create boundary function from timeseries provided in file |
---|
| 191 | #function = file_function(project.boundary_file, domain, verbose=True) |
---|
| 192 | #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) |
---|
| 193 | Br = Reflective_boundary(domain) |
---|
[5095] | 194 | Bd = Dirichlet_boundary([10.,0,0]) |
---|
[5192] | 195 | |
---|
| 196 | |
---|
| 197 | def wave_func(t): |
---|
| 198 | if t < 100: |
---|
| 199 | wave = [(1.2*t/100.0*sin(2*t*pi/10)), 0.0 ,0.0] |
---|
| 200 | else: |
---|
| 201 | wave = [(1.2*sin(2*t*pi/10)), 0.0 ,0.0] |
---|
| 202 | return wave |
---|
| 203 | |
---|
[5095] | 204 | Bwp = Transmissive_Momentum_Set_Stage_boundary(domain, |
---|
[5192] | 205 | wave_func) |
---|
[5154] | 206 | Bwp_velocity = Time_boundary(domain, |
---|
| 207 | lambda t: [(1.2*sin(2*t*pi/10)), |
---|
| 208 | 0.24*sin(2*t*pi/10),0.0]) |
---|
[5095] | 209 | Bws = Transmissive_Momentum_Set_Stage_boundary(domain, |
---|
| 210 | lambda t: [1.2*sin(2*t*pi/10)+1.0*sin(2*t*pi/9.5), |
---|
| 211 | 0.0, 0.0]) |
---|
[5190] | 212 | |
---|
| 213 | if is_structured is True: |
---|
| 214 | domain.set_boundary( {'top': Br, |
---|
| 215 | 'bottom': Br , |
---|
| 216 | 'left' : Br , |
---|
[5192] | 217 | 'right' : Bwp} ) |
---|
[5190] | 218 | else: |
---|
| 219 | domain.set_boundary( {'wall': Br, 'wave': Bwp} ) |
---|
[5093] | 220 | |
---|
[5162] | 221 | |
---|
[5093] | 222 | #------------------------------------------------------------------------- |
---|
| 223 | # Evolve system through time |
---|
| 224 | #------------------------------------------------------------------------- |
---|
| 225 | t0 = time.time() |
---|
| 226 | |
---|
| 227 | for t in domain.evolve(yieldstep, finaltime): |
---|
| 228 | domain.write_time() |
---|
[5192] | 229 | |
---|
[5093] | 230 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
| 231 | print 'finished' |
---|
| 232 | |
---|
[5190] | 233 | points = [[-2,0.0], #-1.8m from SWL |
---|
| 234 | [0,0.0], #0.5m from SWL |
---|
| 235 | [2,0.0], #2m from SWL |
---|
[5093] | 236 | ] |
---|
| 237 | |
---|
| 238 | |
---|
| 239 | #------------------------------------------------------------------------- |
---|
| 240 | # Calculate gauge info |
---|
| 241 | #------------------------------------------------------------------------- |
---|
| 242 | |
---|
[5190] | 243 | if True: |
---|
[5093] | 244 | interpolate_sww2csv(pro_instance.outputdir + basename +".sww", |
---|
| 245 | points, |
---|
[5190] | 246 | pro_instance.outputdir + "depth_dx_"+str(dx)+".csv", |
---|
| 247 | pro_instance.outputdir + "velocity_x_dx_"+str(dx)+".csv", |
---|
| 248 | pro_instance.outputdir + "velocity_y_dx_"+str(dx)+".csv") |
---|
[5093] | 249 | |
---|
| 250 | return pro_instance |
---|
| 251 | |
---|
| 252 | #------------------------------------------------------------- |
---|
| 253 | if __name__ == "__main__": |
---|
[5190] | 254 | main( |
---|
| 255 | is_trial_run=False, |
---|
| 256 | friction=0.0, |
---|
| 257 | is_structured=True, |
---|
[5192] | 258 | outputdir_name='structured_ramp_up2') |
---|