[6160] | 1 | """ |
---|
| 2 | |
---|
| 3 | this test compares the gauge data from a boundary sww file which is precomputed and |
---|
| 4 | stored in this directory, and the same gauge location from a simple model that uses |
---|
| 5 | the same boundary sww file for its boundary conditions |
---|
| 6 | |
---|
| 7 | Basically tests that file_boundary and evolve is working correctly |
---|
| 8 | """ |
---|
| 9 | |
---|
| 10 | #------------------------------------------------------------------------------ |
---|
| 11 | # Import necessary modules |
---|
| 12 | #------------------------------------------------------------------------------ |
---|
| 13 | |
---|
| 14 | # Standard modules |
---|
| 15 | from os import sep,getcwd, access, F_OK, mkdir, getenv |
---|
| 16 | from os.path import dirname, basename,abspath |
---|
| 17 | from shutil import copy |
---|
| 18 | import time, sys, os, tempfile |
---|
| 19 | |
---|
| 20 | # Related major packages |
---|
| 21 | from anuga.shallow_water import Domain,Dirichlet_boundary,File_boundary,Transmissive_boundary, Field_boundary |
---|
| 22 | import Numeric as num |
---|
| 23 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
| 24 | from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files, sww2timeseries, get_data_from_file |
---|
| 25 | # Application specific imports |
---|
| 26 | |
---|
| 27 | #------------------------------------------------------------------------------ |
---|
| 28 | # Copy scripts to time stamped output directory and capture screen |
---|
| 29 | # output to file |
---|
| 30 | #------------------------------------------------------------------------------ |
---|
| 31 | out_dir = dirname(abspath(__file__))+sep |
---|
| 32 | fileName="temp" |
---|
| 33 | |
---|
| 34 | meshes_dir_name = 'small.tsh' # this will be local |
---|
| 35 | |
---|
| 36 | tide = 2.4 |
---|
| 37 | |
---|
| 38 | #-------------------------------------------------------------------------- |
---|
| 39 | # Create the triangular mesh based on overall clipping polygon with a |
---|
| 40 | # tagged |
---|
| 41 | # boundary and interior regions defined in project.py along with |
---|
| 42 | # resolutions (maximal area of per triangle) for each polygon |
---|
| 43 | #-------------------------------------------------------------------------- |
---|
| 44 | |
---|
| 45 | all = [[465184,7764500],[470397,7764510],[470407,7758988],[465195,7758979]] |
---|
| 46 | # N, W, S, E |
---|
| 47 | |
---|
| 48 | create_mesh_from_regions(all, |
---|
| 49 | boundary_tags={'ocean': [ 0],'side': [1, 3], 'back': [2]}, |
---|
| 50 | maximum_triangle_area=50000, |
---|
| 51 | filename=meshes_dir_name, |
---|
| 52 | use_cache=False, |
---|
| 53 | ) |
---|
| 54 | |
---|
| 55 | #------------------------------------------------------------------------- |
---|
| 56 | # Setup computational domain |
---|
| 57 | #------------------------------------------------------------------------- |
---|
| 58 | #print 'Setup computational domain' |
---|
| 59 | domain = Domain( meshes_dir_name, verbose=True) |
---|
| 60 | |
---|
| 61 | from anuga.shallow_water.data_manager import urs2sww |
---|
| 62 | boundaries_dir_name = 'o_test' |
---|
| 63 | |
---|
| 64 | # convert MUX urs files to an SWW file output |
---|
| 65 | urs2sww(boundaries_dir_name,boundaries_dir_name, |
---|
| 66 | mint=9200, maxt= 11200, |
---|
| 67 | fail_on_NaN= False, |
---|
| 68 | verbose=True) |
---|
| 69 | |
---|
| 70 | #------------------------------------------------------------------------- |
---|
| 71 | # Setup initial conditions |
---|
| 72 | #------------------------------------------------------------------------- |
---|
| 73 | |
---|
| 74 | #print 'Start Set quantity' |
---|
| 75 | |
---|
| 76 | domain.set_quantity('elevation', -42.3) |
---|
| 77 | |
---|
| 78 | #print 'Finished Set quantity' |
---|
| 79 | |
---|
| 80 | #------------------------------------------------------ |
---|
| 81 | # Set domain parameters |
---|
| 82 | #------------------------------------------------------ |
---|
| 83 | |
---|
| 84 | domain.set_quantity('stage', tide) |
---|
| 85 | domain.set_quantity('friction', 0.01) |
---|
| 86 | domain.set_name(fileName) |
---|
| 87 | domain.set_datadir(out_dir) |
---|
| 88 | |
---|
| 89 | #------------------------------------------------------------------------- |
---|
| 90 | # Setup boundary conditions |
---|
| 91 | #------------------------------------------------------------------------- |
---|
| 92 | Bf = Field_boundary(out_dir + 'o_test.sww', |
---|
| 93 | domain, mean_stage=tide, use_cache=True, verbose=False) |
---|
| 94 | |
---|
| 95 | #Br = Reflective_boundary(domain) |
---|
| 96 | Bt = Transmissive_boundary(domain) |
---|
| 97 | Bd = Dirichlet_boundary([tide,0,0]) |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | domain.set_boundary({'back': Bt, |
---|
| 101 | 'side': Bd, |
---|
| 102 | 'ocean': Bf |
---|
| 103 | }) |
---|
| 104 | #---------------------------------------------------------------------------- |
---|
| 105 | # Evolve system through time |
---|
| 106 | #---------------------------------------------------------------------------- |
---|
| 107 | |
---|
| 108 | t0 = time.time() |
---|
| 109 | |
---|
| 110 | for t in domain.evolve(yieldstep = 60, finaltime = 1920): |
---|
| 111 | domain.write_time() |
---|
| 112 | domain.write_boundary_statistics() |
---|
| 113 | |
---|
| 114 | |
---|
| 115 | #---------------------------------------------------------------------------- |
---|
| 116 | #Gets timeseries from boundary sww and evolved sww |
---|
| 117 | #---------------------------------------------------------------------------- |
---|
| 118 | home = getenv('INUNDATIONHOME') #Sandpit's parent dir |
---|
| 119 | data = 'data' |
---|
| 120 | state = 'western_australia' |
---|
| 121 | scenario_name = 'dampier.sww' |
---|
| 122 | |
---|
| 123 | scenario = 'dampier_tsunami_scenario_2006' |
---|
| 124 | an = 'anuga' |
---|
| 125 | bo = 'boundaries' |
---|
| 126 | |
---|
| 127 | run_time = 'blank' |
---|
| 128 | production_dirs = {run_time: 'URS evolved data'#, |
---|
| 129 | } |
---|
| 130 | |
---|
| 131 | topo = 'topographies' |
---|
| 132 | out = 'outputs' |
---|
| 133 | urs = 'urs' |
---|
| 134 | gridded = '1_10000' |
---|
| 135 | |
---|
| 136 | |
---|
| 137 | gauge_boundary_filename = 'gauges_time_series_b_near_top.csv' |
---|
| 138 | gauge_evolved_filename = 'gauges_time_series_near_top.csv' |
---|
| 139 | |
---|
| 140 | boundary_dir_filename = os.path.join(out_dir,gauge_boundary_filename) |
---|
| 141 | #print'boundary_dir_filename',boundary_dir_filename |
---|
| 142 | |
---|
| 143 | evolved_dir_filename= os.path.join(out_dir,gauge_evolved_filename) |
---|
| 144 | |
---|
| 145 | swwfiles = {} |
---|
| 146 | swwfile = out_dir + fileName + '.sww' |
---|
| 147 | swwfiles[swwfile] = run_time |
---|
| 148 | #---------------------------------------------------------------------------- |
---|
| 149 | #get timeseries data from evolved sww file |
---|
| 150 | #---------------------------------------------------------------------------- |
---|
| 151 | texname, elev_output = sww2timeseries(swwfiles, |
---|
| 152 | out_dir+sep+"gauges.txt", |
---|
| 153 | # out_dir+sep+"gauges.csv", |
---|
| 154 | production_dirs, |
---|
| 155 | report = False, |
---|
| 156 | plot_quantity = ['stage', 'xmomentum', 'ymomentum'], |
---|
| 157 | surface = False, |
---|
| 158 | time_min = None, |
---|
| 159 | time_max = None, |
---|
| 160 | title_on = False, |
---|
| 161 | use_cache = True, |
---|
| 162 | verbose = False) |
---|
| 163 | |
---|
| 164 | swwfiles = {} |
---|
| 165 | swwfile = out_dir + boundaries_dir_name + '.sww' |
---|
| 166 | swwfiles[swwfile] = run_time |
---|
| 167 | #print"swwfiles",swwfiles,"shallow_water" |
---|
| 168 | |
---|
| 169 | #---------------------------------------------------------------------------- |
---|
| 170 | #get timeseries data from boundary sww file |
---|
| 171 | #---------------------------------------------------------------------------- |
---|
| 172 | texname, elev_output = sww2timeseries(swwfiles, |
---|
| 173 | out_dir+sep+"boundary_gauges.txt", |
---|
| 174 | # out_dir+sep+"boundary_gauges.csv", |
---|
| 175 | production_dirs, |
---|
| 176 | report = False, |
---|
| 177 | plot_quantity = ['stage', 'xmomentum', 'ymomentum'], |
---|
| 178 | surface = False, |
---|
| 179 | time_min = None, |
---|
| 180 | time_max = None, |
---|
| 181 | title_on = False, |
---|
| 182 | use_cache = True, |
---|
| 183 | verbose = False) |
---|
| 184 | |
---|
| 185 | #---------------------------------------------------------------------------- |
---|
| 186 | #get_data_from file returns a e_data which is an array containing the |
---|
| 187 | #data from the file it read |
---|
| 188 | #---------------------------------------------------------------------------- |
---|
| 189 | |
---|
| 190 | e_header, e_data = get_data_from_file(evolved_dir_filename) |
---|
| 191 | |
---|
| 192 | # assign columns from array to single vector(arrays) |
---|
| 193 | e_time = e_data[:,0] |
---|
| 194 | e_stage = e_data[:,1] |
---|
| 195 | e_momentum = e_data[:,2] |
---|
| 196 | e_speed = e_data[:,3] |
---|
| 197 | e_elevation = e_data[:,4] |
---|
| 198 | |
---|
| 199 | #---------------------------------------------------------------------------- |
---|
| 200 | #get_data_from file returns a b_data which is an array containing the |
---|
| 201 | #data from the file it read |
---|
| 202 | #---------------------------------------------------------------------------- |
---|
| 203 | |
---|
| 204 | b_header, b_data = get_data_from_file(boundary_dir_filename) |
---|
| 205 | |
---|
| 206 | # assign columns from array to single vector(arrays) |
---|
| 207 | b_time = b_data[:,0] |
---|
| 208 | b_stage = b_data[:,1] |
---|
| 209 | b_momentum = b_data[:,2] |
---|
| 210 | b_speed = b_data[:,3] |
---|
| 211 | b_elevation = b_data[:,4] |
---|
| 212 | |
---|
| 213 | #---------------------------------------------------------------------------- |
---|
| 214 | #compares the 2 models, make the vector of data from the boundary sww |
---|
| 215 | #file have the same time interval as the evolved sww vector. this |
---|
| 216 | #allows them to be compared with the allclose statment at the end |
---|
| 217 | #---------------------------------------------------------------------------- |
---|
| 218 | |
---|
| 219 | j=0 |
---|
| 220 | k=0 |
---|
| 221 | #from boundary sww |
---|
| 222 | b_sample = [] |
---|
| 223 | #from evolved sww |
---|
| 224 | e_sample = [] |
---|
| 225 | for i in range(len(b_time)): |
---|
| 226 | # if j<(len(e_time)) and b_time[i] == e_time[j]: |
---|
| 227 | if j<(len(e_time)) and k<(len(e_time)) and b_time[i] == e_time[j]: |
---|
| 228 | b_sample.append(float(tide+b_stage[i])) |
---|
| 229 | e_sample.append(float(e_stage[k])) |
---|
| 230 | # if k <len(e_time): print 'time e equal b:', b_time[i],i, j, b_stage[i], b_sample[i], e_stage[k],(len(e_time)-1) |
---|
| 231 | j = j +1 |
---|
| 232 | k = k +1 |
---|
| 233 | |
---|
| 234 | |
---|
| 235 | e_stage = e_stage.tolist() |
---|
| 236 | e_stage.pop() |
---|
| 237 | |
---|
| 238 | |
---|
| 239 | print "test successful" |
---|
| 240 | print 'fileName',fileName+'.sww' |
---|
| 241 | os.remove(meshes_dir_name) |
---|
| 242 | #---------------------------------------------------------------------------- |
---|
| 243 | # Compares the two time series |
---|
| 244 | #---------------------------------------------------------------------------- |
---|
| 245 | assert num.allclose (b_sample, e_sample, 0.5, 0.5) |
---|
| 246 | |
---|
| 247 | |
---|
| 248 | |
---|
| 249 | |
---|
| 250 | |
---|
| 251 | |
---|
| 252 | |
---|
| 253 | |
---|
| 254 | |
---|
| 255 | |
---|
| 256 | |
---|
| 257 | |
---|
| 258 | |
---|
| 259 | |
---|
| 260 | |
---|